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Abstract. Neutrino transport and neutrino interactions in dense matter play a crucial role in stellar core collapse, 
supernova explosions and neutron star formation. Here we present a detailed description of a new numerical code 
for treating the time and energy dependent neutrino transport in hydrodynamical simulations of such events. 
The code is based on a variable Eddington factor method to deal with the integro-differential character of the 
Boltzmann equation. The moments of the neutrino distribution function and the energy and lepton number ex- 
change with the stellar medium are determined by iteratively solving the zeroth and first order moment equations 
in combination with a model Boltzmann equation. The latter is discretized on a grid of tangent rays. The in- 
tegration of the transport equations and the neutrino source terms is performed in a time-implicit way. In the 
present version of the program, the transport part is coupled to an explicit hydrodynamics code which follows the 
evolution of the stellar plasma by a finite-volume method with piecewise parabolic interpolation, using a Riemann 
solver for calculating the hydrodynamic states. The neutrino source terms are implemented in an operator-split 
step. Neutrino transport and hydrodynamics can be calculated with different spatial grids and different time 
steppings. The structure of the described code is modular and offers a high degree of flexibility for an application 
to relativistic and multi-dimensional problems at different levels of refinement and accuracy. We critically evaluate 
results for a number of test cases, including neutrino transport in rapidly moving stellar media and approximate 
relativistic core collapse, and suggest a path for generalizing the code to be used in multi-dimensional simulations 
of convection in neutron stars and supernovae. 
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1. Introduction 

■ When the core of a massive star collapses to a neutron 
, star, a huge amount of gravitational binding energy is re- 
leased mainly in the form of neutrinos which are abun- 
dantly produced by particle reactions in the hot and dense 
plasmaj. In facr, rhe emission of neurrinos determines the 



sequence of dynamical events which precede the death of 
the star in a supernova explosion. Electron neutrinos from 
electron captures on protons and nuclei reduce the elec- 
tron fraction and the pressure and thus accelerate the con- 
traction of the stellar iron core to a catastrophic implo- 
sion. The position of the formation of the supernova shock 
at the moment of core bounce depends on the degree of 
deleptonization during the collapse phase. The outward 
propagation of the prompt shock is severely damped by 
energy losses due to the photodisintegration of iron nu- 
clei, and finally stopped only a few milliseconds later when 



tron neutrinos. These neutrinos are created in the shock- 
heated matter and leave the star as soon as their diffusion 
is faster than the expansion of the shock. "Prompt" su- 
pernova explosions generated by a direct propagation of 
the hydrodynamical shock have been obtained in simula- 
tions only when the stellar iron core is very small (Baron & 
Cooperstcin 199(\ ) and/or the equation of state of neutron- 



rich nuclear matter is extraordinarily soft (Baron et al 



1985, 1987) 



At a later stage (roughly a hundred milliseconds after 
the shock formation) the situation has changed. The tem- 
perature behind the stalled shock has dropped such that 
increasingly energetic neutrinos diffusing out from deeper 
layers start to transfer energy to the stellar gas around 
the nascent neutron star. If this energy deposition is suf- 
ficiently strong, the stalled shock can be revived and a 



"delayed" explosion can be triggered (Wilson 1985; Bethc 



additicnal energy is lost in a luminous outburst of elec- & Wilson 1985; the idea that neutrinos provide the energy 
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of the supernova explosion was originally brought up by 
polgate fc White 1966). A small fraction of less than one 
per cent of the energy released in neutrinos can account for 
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the kinetic energy of a typical type II supernova. This ex- & Mayle 1993) develops on a time scale of tens of millisec 



planation is currently favored for the explosion of massive 
stars in the mass range between about 10 Mq and roughly 
25 Mq. It is supported by numerical simulations and ana- 
lytic considerations. A finally convincing hydrodynamical 
simulation, however, is still missing (a summary of our 
current understanding of the explosion mechanism can be 



onds after core bounce (Keil et al. 1996; Keil 1997; also 



lanka et al. 2001). The models were constrained to a sim- 
ulation of the neutrino cooling of the nascent neutron star, 
but it must be expected that the convective enhancement 
of the neutrino luminosity, which becomes appreciable af- 



found in lanka 2001). The measurement of a neutrino sig- 
nal from a Galactic supernova would offer the most direct 
observational test for our theoretical perception of the ou- 



ter about 200 milliseconds (see also Janka et al. 2001), 
can have important consequences for the revival of the 
stalled supernova shock ( [Burrow^ 1987). Ledoux unsta- 
ble regions were also detected in spherical cooling models 



set of the explosion. Due to the central role of neutrinos of newly formed neutron stars ( 


Burrows 1987 


; I 


3urrows 


in supc rnovae, the neutrino transport deserves particular & Lattimer 


1988; 


Pons et al. 


1999; 


Miralles et a" 


. 2000). 



attention in numerical models. 



1.1. Brief review of the status of supernova models 

Current hydrodynamic models of neutrino-driven super- 
nova explosions leave an ambiguous impression and have 
caused confusion about the status of the field outside the 
small community of supernova modelers. Simulations by 
different groups seem to be contradictory because some 
models show successful explosions by the neutrino-driven 
mechanism whereas others have found the mechanism to 
fail. 

Wilson and collaborators have performed successful 



simulations for more than 15 years now (e.g.. 


Wilsoi-^ 


1985; 


Wilson et al 


1986; 


Wilson & Maylc 


1988, 


1993; 1 


Mayle 



et al. 1993 



Nevertheless, their significance is controversial, because 
Bruenn et al. ( 1995) in sp herically symmetric models and 
Mezzacappa et al. ( 1998a ) in two-dimensional simulations 
observed convection inside the neutrinosphere only as a 
transient phenomenon during a relatively short period af- 
ter core bounce. 

The question is therefore undecided whether convec- 
tion plays an important role in newly formed neutron 
stars and if so, what its implications for the super- 
nova explosion are. Unfortunately, the various calcula- 
tions were performed with different treatments of the neu- 
trino transport, one-dimensional vs. two-dimensional hy- 
drodynamics, general relativistic gravitational potential 
vs. Newtonian gravity, or even an inner boundary condi- 
tion to replace the central, dense part of the neutron star 



Totani et al. 1998) Their models were cal- 
culated in spherical symmetry, but shock revival was ob- 
tained by making the assumption that the neutrino lu- 
minosity is boosted by neutron-finger instabilities in the 
nascent neutron star (Wilson & Mayk 198^ , 1993). In 



core (Mezzacappa et al. 



however, where Keil et al. 



1998a). It is this inner region, 



(1996) found convection to de- 
velop roughly 100 milliseconds after bounce. Convection 
at larger radii and thus closer below the neutrinosphere 
had indeed died out within only 20-30 milliseconds after 



shock formation in agreement with the results of Bruenn 



to typically observed values (of the order of 10^^ erg) re 
quire pions to be abundant in the nuclear matter already 
at moderately high densities ( pVIayle et al. 1993). The cor- 
responding equation of state (EoS) with pions leads to 
higher core temperatures and the emission of more ener- 
getic neutrinos. Both the convectively boosted luminosi- 
ties and the harder spectra enhance the neutrino heat- 
ing behind the stalled shock. The relevance of neutron- 
finger instabilities for efficient energy transport on large 
scales, however, has not been demonstrated by direct sim- 
ulations. The existence of neutron-finger instabilities (in 



fact, ii| then- calculations explosion energies comparable ^t al.| (p^95|). Future studies with a better and more con- 



sistent handling of the different aspects of the physics are 
definitely needed to clarify the situation. 

A second hydrodynamically unstable region develops 
exterior to the nascent neutron star in the neutrino-heated 
layer behind the stalled supernova shock. Convective over- 



in cases which otherwise fail ( 


Herant et al. 


1992, 


1994; 


Janka & Miiller 


1995, 


1996; 


Burrows et al. 


1990. In two- 



dimensional supernova models computed recently by Fryei 
, [Fryer fc Heg"a] (|2GG0| ), and [Fryer et al.| (|00|) these 



the sense of the definition introduced by Wilson fc Mayle 



1995) is indeed questioned by other investigations (Bruenn 



fc Dineva 1996 ) and might be a consequence of specific 
properties of the high-density EoS or the treatment of 
the neutrino transport by Wilson and collaborators. Also 
the implementation of a pionic component in their EoS 
is at least controversial and not in agreement with other, 
more conventional descriptions of nuclear matter (Pethick 
& Pandharipande, personal communication). 

Rather than finding neutron-finger instabilities, two- 
dimensional hydrodynamic simulations have shown that 
regions inside the nascent neutron star can exist where 



Ledoux or quasi-Ledoux convection (as defined by Wilson 



( [199' 

hydrodynamical instabilities in the postshock region are 
crucial for the success of the neutrino-driven mechanism, 
because they help transporting hot gas from the neutrino- 
heating region directly to the shock, while downftows si- 
multaneously carry cold, accreted matter to the layer of 
strongest neutrino heating where a part of this gas readily 
absorbs more energy from the neutrinos. The existence of 
this multi-dimensional phenomenon seems to be generic 
for the situation which builds up in the stellar core some 
time after shock formation. It is therefore no matter of 
dispute. 

All simulations showing explosions as a consequence 
of post-shock convection, however, have so far been per- 
formed with a strongly simplified treatment of the cru- 
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cial neutrino physics, e.g., with grey, flux-limited diffu- 
sion which is matched to a "light bulb" description at 
some "low" value of the optical depth (Herant et al. 1994; 
Burrows et al. 1995). Alternatively, an inner boundary 



near the neutrinosphere had been used where the spec- 
tra and luminosities of neutrinos were prescribed to pa- 
rameterize our ignorance and the potential uncertainties 
of the exact properties of the neutrino emission from the 
newly formed neutron star (Janka & Miiller 1995, 1996| ). 



It is therefore not clear whether the instabilities and their 
associated effects in fully self-consistent and more ac- 
curate simulations will be strong enough to cause suc- 
cessful explosions. Doubt s in that respect were raised by 
Mezzacappa et al. ( 1998b ), whose two-dimensional models 
showed convective overturn in the neutrino-heating region 
but still no explosion. Mezzacappa et al. ( |l998b ) combined 
two-dimensional hydrodynamics with neutrino transport 
results obtained by multi-group flux-limited diffusion in 
spherical symmetry. Although not self-consistent, their ap- 
proach is nevertheless an improvement compared to pre- 
vious treatments of the neutrino physics by other groups. 

All computed models bear some pieces of truth. 
Essentially they demonstrate the remarkable sensitivity of 
the supernova dynamics to the different physical aspects of 
the problem, in particular the treatment of neutrino trans- 
port and neutrino-matter interactions, the properties of 
the nuclear EoS, multi-dimensional hydrodynamical pro- 
cesses, and general relativity. Considering the huge energy 
reservoir carried away by neutrinos, the neutrino-driven 
mechanism appears rather inefficient (the often quoted 
value of 1% efliciency, however, misjudges the true situa- 
tion, because neutrinos can transfer between 5% and 10% 
of their energy to the stellar gas during the critical period 
of shock revival; see Janka 2001). Nevertheless, neutrino- 



driven explosions are "marginal" in the sense that the 
energy of a standard supernova explosion is of the same 
order as the gravitational binding energy of the ejected 
progenitor mass. The final success of the supernova shock 
is the result of different physical processes which compete 
against each other. On the one hand, neutrino heating 
tries to drive the shock expansion, on the other hand en- 
ergy losses, e.g., by neutrinos that are reemitted from the 
inward flow of neutrino-heated matter which enters the 
neutrino-cooling zone below the heating layer, extract en- 
ergy and thus damp the shock revival. It is therefore not 
astonishing that different approximations or descriptions 
for one or more physical components of the problem can 
decide between an explosion or failure of a simulation. 

None of the current supernova models includes all rel- 
evant aspects to a satisfactory level of accuracy, but all 
of these models are deflcient in one or more respects. 
Wilson and collaborators obtain explosions, but their in- 
put physics is unique and cannot be considered as gen- 
erally accepted 



Two-dimensional (Herant 
Burrows etH] |1995| ; |Fryei) |1999| ; [Fryer 



et al 



Heger 



1994; 



2000; 



Fryer et al. 2002) and three-dimensional (Fryer, personal 
communication) models show explosions due to strong 
postshock convection, but the neutrino transport and 



neutrino-matter interactions are handled at a level of ac- 
curacy which falls back behind the most elaborate treat- 
ments that have been applied in spherical symmetry. The 
sensitivity of the outcome of numerical calculations to de- 
tails of the neutrino transport demands improvements. 
Self-consistent multi-dimensional simulations with a so- 
phisticated and quantitatively reliable treatment of the 
neutrino physics have yet to be performed. 

1.2. Boltzmann solvers for neutrino transport 



ativistic ( Liebcndorfer et al. 2001b) hydrodynamical sim- 



Most recently, spherically symmetric Newtonian (Rampp 
|fc Janka 2000; Mezzacappa et al. 2001), and general rel 



ulations of stellar corc-coUapse and post-bounce evolu- 
tion including a Boltzmann solver for the neutrino trans- 
port have become possible. Although the models do not 
yield explosions, they must be considered as a major 
achievement for the modeling of supernovae. Before these 
calculations only the collapse phase of the stellar core 
had b een investigated with so l ving t he Boltzmann equa- 



tion ( Mezzacappa fc Bruenn 1993c ). Boltzmann trans- 



port has now superseded multi-group flux-limited diffu- 



sion ( 


Arnctt 


1977; Bowers & Wilson 


1982 


Bruenn 


1985 


Myra et al. 


1987; 


Baron et al. 


198£: 


Bruenn 


1989a 


,b 


1993 


Cooperstcin & Baron 


1992 


) as the most elaborate treat- 



mcnt of neutrinos in supernova models. 

The differences between both methods in dynamical 
calculations have still to be figured out in detail, but 
the possibility of solving the Boltzmann equation removes 
imponderabilities and inaccuracies associated with the 
use of fiux-limiters in particular in the region of semi- 
transparency, where neutrinos decouple from the stellar 
background and also deposit the energy for an explo- 



sion 


( Janka 


1991, 


1992; 


Messer et al. 


1998; 


Yamada et al. 


1999 


). For the first time the transport can now be han- 



died at a level of sophistication where the technical treat- 
ment is more accurate than our standard description of 
neutrino-matter interactions, which includes various ap- 
proximations and simplifications (for an overview of the 
status of the ha ndling of neutrino-nucleon interactions, see 



Horowitz 2002 and the references therein). 

In this paper we describe our new numerical code for 
solving the energy and time dependent Boltzmann trans- 
port equation for neutrinos coupled to the hydrodynamics 
of the stellar medium. First results from supernova calcu- 
lations with this code were published before (Rampp & 
panka 200C| ) , but a detailed documentation of the method 
will be given here. It is based on a variable Eddington fac- 
tor technique where the coupled set of Boltzmann equa- 
tion and neutrino energy and momentum equations is iter- 
ated to convergence. Variable Eddington moments of the 
neutrino phase space distribution arc used for closing the 
moment equations, and the integro-diffcrcntial character 
of the Boltzmann equation is tamed by using the zeroth 
and first order angular moments (neutrino density and 
flux) in the source terms on the right hand side of the 
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Boltzmann equation. This numerical approach is funda- 



mentally different from the so-called Sn methods (( 


I!arlson 


1967 


Yueh & Buchler 


1977; 


Mezzacappa & Bruenn 


1993a; 



Yamada et al. 1999) which employ a direct discretization 
of the Boltzmann equation in all variables including the 
dependence on the angular direction of the radiation prop- 
agation. 

Some basic characteristics of our code are similar to el- 
ements described by Burrows et al. ( [2000 ) . Different from 



the latter paper we shall discuss the details of the nu- 
merical implementation of the transport scheme and its 
coupling to a hydrodynamics code (Elampp 2000| ), which 
in our case is the PROMETHEUS code with the po- 
tential of performing multi-dimensional simulations. The 
variable Eddington factor technique was our method of 
choice for the neutrino transport because of its modularity 
and flexibility, which offer significant advantages for a gen- 
eralization to multi-dimensional problems. We shall sug- 
gest and motivate corresponding approximations which we 
consider as reasonable in the supernova case, at least as a 
first step towards multi-dimensional hydrodynamics with 
a Boltzmann treatment of the neutrino transport. 

This paper is arranged as follows: In Sect. ^ the equa- 
tions of radiation hydrodynamics are introduced. Sect. ^ 
provides a general overview of the numerical methods used 
to solve these equations and contains details of their prac- 
tical implementation. Results for a number of idealized 
test problems as well as applications of the new method 
to the supernova problem are presented in Sect. ^ A sum- 
mary will be given in Sect ^. In the Appendix the numeri- 
cal implementation of neutrino opacities and the equation 
of state used for core-collapse and supernova simulations 
is described. 



2. Radiation-hydrodynamics — basic equations 

2.1. Hydrodynamics and equation of state 

For an ideal fluid characterized by the mass density p, the 
Cartesian components of the velocity vector (wi, W2, wa)"^, 
the specific energy density e = e+^v'^ and the gas pressure 
p, the Eulerian, nonrelativistic equations of hydrodynam- 
ics in Cartesian coordinates read (sum over i implied): 

dtp + d^ipv,) = , (1) 
dtipvk) + d,{pv,vk + d,kP) = -p^fe^^^^'V Quk , (2) 
dtipe) + d,{{pe+p}vi) = Qe+ v,Qm, , (3) 

where (j)Ncwt (j^gnotes the Newtonian gravitational poten- 
tial of the fluid, which can be determined by the Poisson 
equation 9^5'$^™* = inG p (G is Newton's constant). 
<5m and Qe are the neutrino source terms for momen- 
tum transfer and energy exchange, respectively, S^f^ is the 
Kronecker symbol, and di := d /dx^ is an abbreviation for 
the partial derivative with respect to the coordinate . In 
order to describe the evolution of the chemical composi- 
tion of the (electrically neutral) fluid, the hydrodynamical 
equations are supplemented by a conservation equation for 



the electron fraction 1^, 

dt{pY,) + d,{pY,Vi) 



(4) 



where the source term Qn describes the change of the net 
electron number density (i.e. the density of electrons mi- 
nus that of positrons) due to emission and absorption of 
electron-flavour neutrinos. Unless nuclear statistical equi- 
librium (NSE) holds, an equation like (^ has to be solved 
for the mass fraction Xk of each of the A^nuc individual (nu- 
clear) species. In NSE, Xk — Xk{p,T,Yc) is determined 
by the Saha equations. 

An equation of state is invoked in order to express 
the pressure as a function of the independent thermody- 
namical variables, i.e., p = p{p,T,Yc), if NSE holds, or 
p = p{p, T, Yc, {Xk}k=i jY ) otherwise (see Appendix ^ 
for the numerical handling of the equation of state). 

In the following we will employ spherical coordinates 
and, unless otherwise stated, assume spherical symmetry. 

2.2. Equations for the neutrino transport 

2.2.1. General relativistic transport equation 

Lmdquist ( |1966 ) derived a covariant transfer equation and 
specialized it for particles of zero rest mass interacting in 
a spherically symmetric medium supplemented with the 
comoving frame metric (a is a Lagrangian coordinate) 
= _e2*(*.«)c2dt2 + e2A(*.«)da2 + R(t, a)W 
The "Lindquist-equation" , which describes the evolu- 
tion of the specific intensity I as measured in the comoving 
frame of reference, reads: 

1 I - u? d 

-Vtl + pVal + T I 
c R op 



d_ 
dp 

d_ 
'd'e 



,1 



( 



u 



VtA + pVa^]l 



(3 - A^');^ + (1 + m')-AA + 2pVa<P ]I = G, 



(5) 



where we use the classical abbreviations r{t, a) := VaR 
and U(t,a) := c'^VtR with Va := e~'^da and Vt := 
e~*9(. The latter definition has been used in Eq. (||) also 
in order to emphasize that the time derivative has to be 
taken at fixed Lagrangian coordinate a. 

The functional dependences of the metric functions 
^(f, a), A(t, a), R{t,a), the specific intensity X{t,a,e, p), 
and the collision integral G{t, a, e, p) were suppressed for 
brevity. Momentum space is described by the coordinates 
e and p, which are the energy and the cosine of the angle 
of propagation of the neutrino with respect to the radial 
direction, both measured in the locally comoving frame 
of reference. Note that the opacity x and the emissivity 
77, and thus the collision integral C = rj — yX in general 
depend also explicitly on momentum-space integrals of I, 
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which makes the transfer equation an integro-partial dif- 
ferential equation. Examples of the actual computation of 
the collision integral for a number of interaction processes 
of neutrinos with matter can be found in Appendix 

2.2.2. 0{v/c) transport equations 

In general, the metric functions <f>(t, a), A(i, a) and R{t, a) 
have to be computed numerically from the Einstein field 
equations. When working to order 0{f3:=v/c) and in a 
flat spacetime (usually called the "Newtonian approxima- 
tion"), it is however possible to express these functions 
analytically in terms of only the velocity field and its 
first time derivative (the fluid acceleration). Details of the 
derivation can be found in Castoi| (1972). Alternatively 
one can simply reduce the special relativistic transfer 
equation ( |Mihalas| [T980| ) to order 0{v/c). 

This transfer equation, together with its angular mo- 
ment equations of zeroth and first order reads (e.g., 
Mihalas fc Mihalas||l984 see also|Lowrie et al||200l[): 



,ld_ 

'cdt 



4)1- 

or 
d_ 

d_ 



or 



r dfj, 



s((i-^^)^+M 



[3 ,d(3 ■ 



+ ((3-^2)^ + (l + M^) 



dr 
dr 



2 9/3\^ 



C\ (6) 



(-- 



d_ 

We 



r or 



f^{3J-K) + ^{J + K) 
r or 



cm" 



,\d_ 

'cdt 



d 



„ \ d , K-J 
r or c ot 



dr 
d_ 



21? + - 

or r 



c dt 



(8) 



For reference we also write down the transformations 
(correct to 0{P)) which allow one to relate the frequency- 
integrated moments in the comoving ( "Lagrangian" ) and 
in the inertial ("Eulerian") frame of reference (indicated 
by the superscript "Eul" ) . 

Ji="' = J+2/3i7, 



= H+13 (J + K), 



(10) 



K 



Eul 



K+2(3H. 



Eq. (|10|) can easily be deduced from a Lorentz- 
transformation of the radiation stress-energy tensor (e.g., 
Mihalas & Mihalas 1984). In principle also relations for 
monochromatic moments can be derived by transforming 
the specific intensity X, the angle cosine /i and the energy 
e, which, however, leads to more complicated expressions. 

The system of Eqs. (^-||) is coupled to the evolution 
equations of the fluid (Eqs. in spherical coordinates 
and symmetry) by virtue of the definitions of the source 
terms 



Qe = -47r 
— 



deC(°)(e), 
deC«(e), 
47rmB / deC(")(e) 



(11) 
(12) 
(13) 



where me denotes the baryonic mass, and C(°)(e) 
e-iC(o)(e). 



Recently, Lowrie et al. ( 2001 ) emphasized the funda- 
mental significance of a term 



p,-'4i 

c dt 



(14) 



on the left hand side of the special relativistic transport 
equation. This term has traditionally been dropped in de- 
riving the C'(t;/c)-approximation (Eq. ^) by assuming the 
time variation of all dependent variables (like, e.g., T) to 
be on a fluid time scale (given by l/v, where I is the char- 
acteristic length scale of the system and v a typical fluid 
velocity). In this case the ter m given by Eq. (|l^ is of o rder 
v'^/c^ and can be omitted ( Mihalas fc Mihalas 1984). If, 
on the other hand, appreciable temporal variations of, e.g., 
the specific intensity occur on radiation time scales (given 
by l/c), E q. (^) is no longer correct to 0{v/c) (Lowrie 



where, lin spherical symmetry, the angular moments of the 
specific intensity are given by 

+1 

{J,H,K,L,... }{t, r, e) := ^ / dA* fi^"'^'^'- ^I{t, r, e, m) • 

(9) 

Note that in Eqs. (^-^ all physical quantities, in par- 
ticular also the collision integral and its angular moments 
C'^'^-' {t, r, e) := ^ J^^ d/i ^^C{t, r, e, fi) are measured in the 
comoving frame, but the choice of coordinates (r, t) is 



et al. 2001) 



In core-collapse supernova simulations carried out so 
far, the dynamics of the stellar fluid presumably was not 
affected by neglecting the term in Eq. (^. However, our 
tests with Eq. (||) , including the additional time derivative 
of Eq. (^) and the corresponding changes in the moment 
equations (Eqs. |7| ^), have shown that the neutrino signal 
computed in a supernova simulation is indeed altered com- 
pared to the traditional treatment. Wc will therefore take 
the term of Eq. ( p^ ) into account in future simulations. 

For calculating nonrelativistic problems Mihalas & 



Eulerifn. The simple replacement d/dt + vd/dr D/Di 



yields the conversion to Lagrangian coordinates. 



Mihalas (1984) suggested a form of the radiation mo- 
mentum equation (Eq. ||) , in which all velocity-dependent 
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terms except for the (3d / dr-ierra in the first hne of Eq. (||) 
are dropped. When the velocities become sizeable, how- 
ever, it may be advisable to solve the momentum equation 
in its general form (E q. p|). Doing so, we in deed found that 
the terms omitted by Mihalas fc Mihalas] ( |1984| ) can have 



an effect on the solution of the neutrino transport in su- 
pernovae, in particular on the neutrino energy spectrum. 



3. Numerical implementation 

3.1. Overview of the basic numerical strategy 

For solving the neutrino radiation hydrodynamics problem 
we employ the well-known "operator splitting" -approach, 
a popular method to make large problems computation- 
ally tractable by considering them as a combination of 
smaller subproblems. The coupled set of equations of radi- 
ation hydrodynamics (Eqs. ^-^) is split into a "hydro- 
dynamics part" and a "neutrino part" , which are solved 
independently in subsequent ( "fr actional" ) steps. In the 
"hydrodynamics part" (Sect. 3^) a solution of the equa- 



tions of hydrodynamics including the effects of gravity is 
computed without taking into account neutrino-matter in- 
teractions. 

In what we call the "neutrino part" (Sects. 3.3-3.4) 



the coupled system of equations describing the neutrino 
transport and the changes of the specific internal energy 
e and the electron fraction Yc of the stellar fluid due to 
the emission and absorption of neutrinos are solved. The 
latter changes are given by the equations 



S An 
St p Ja 



S ^ Airm 
St " P Jo 



i.e{i.e,Pe,...} 

de (c^^He) ~ Ci:\e)) . (16) 



Here, the notation SJSt indicates that the full problem as 
given, e.g., by Eq. (jj), is split into the equations dt{pYc) + 
r~'^dr {r^pYc v) = and 5t Yc = Qn/p. 

Within the neutrino part, we calculate in a first step 
the transport of energy and electron-type lepton number 
by neutrinos and the corresponding exchange with the 
stellar fluid for electron neutrinos {vc) and antineutrinos 
(Dc) only, which also means that the sum in Eq. ( |l5| ) is 
restricted to u € {vq^Dc}. The /i and r neutrinos and their 
antiparticles, all represented by a single species ("i^x") are 
handled together with the remaining terms of the sum in 
Eq. ( p^ ) in a second step. 

For the neutrino transport we use the so-called "vari- 
able Eddington factor" method, which determines the 
neutrino phase-space distribution by iteratively solving 
the radiation moment equations for neutrino energy and 
momentum coupled to the Boltzmann transport equa- 
tion. Closure of the set of moment equations is achieved 
by a variable Eddington factor calculated from a for- 
mal solution of the Boltzmann equation, and the integro- 
differential character of the latter is tamed by making use 



of the integral moments of the neutrino distribution as ob- 
tained from the moment equations. The method is similar 
to the one described by Burrows ct al.| ( |2000| ). 



3.2. The hydrodynamics code 

For the integration of the equations of hydrodynamics 
we employ the Newtonian finite-vo lume hydrodynamic s 
code PROMETHEUS developed by [Fryxell eFall ( |l989|) . 
PROMETHEUS is a direct Eulerian, time-explicit im- 
ple mentation of the Piecew i se Pa rabolic Method (PPM) 
of Colella fc Woodward ( 1984 ), which is a second- 
order Godunov scheme employing a Riemann solver. 
PROMETHEUS is particularly well suited for following 
discontinuities in the fiuid flow like shocks or boundaries 
between layers of different chemical composition. It is ca- 
pable of tackling multi-dimensional problems with high 
computational efficiency and numerical accuracy. 

The version of PROMETHEUS used in our radiation 



hydrodynamics code was provided to us by |K!eil| (|1997|) . It 
offers an optional generalization of the Newtonian gravita- 
tional potential to include general relativistic corrections. 
The hydrodynamics code was considerably updated and 
augmented by K. Kifonidis who added a simplified version 
of the "Consistent Mul tifluid Advection (CMA)" method 
( Plcwa fc Miiller 1999) which ensures an accurate advec- 
tion of the individual chemical components of the fluid. 
In order to avoid spurious oscillations arising in multi- 
dimensional simulations, when strong shocks are aligned 
with one of the coordinate lines (the so-called "odd-even 
decoupling" phenomenon), a hybrid scheme was imple- 
mented dKifonidi^ |2000| ; [Plewa fc Miillci] ^00 1| ) which, in 



the vicinity of such shocks, selectively switches from the 
original PPM method to the more diffusive HLLE solver 



of Einfeldt (1988 



3.3. Multi-group formulation of the 0{v/c) moment 
equations 

In what follows we describe the numerical implementa- 
tion of the "Newtonian" limit of the transport equations 
including 0(v/c)-terms (cf. Sect. 2.2.2). The general rela- 
tivistic version of the code is based on almost exactly the 
same considerations. 



3.3.1. General considerations 

In order to construct a conservative numerical scheme for 
the system of moment equations, we integrate the moment 
equations (Eqs. 0, ||) over a zone of the numerical (r, e)- 
mesh. 



Spatial discretization: After having performed the integral 
over the volumes Ay of the radial mesh zones, the mo- 
ment equations can be rewritten as evolution equations 
for the volume-integrated moments (e.g., ^^yJdV). In 
case of a moving radial grid, where the volume AV^ of a 
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grid cell is time dependent, one has to apply the "moving 
grid transport theorem" (see e.g., [Miillei 1991) in order to 
interchange the operators Dj and J^y dV. For the special 
case of a Lagrangian grid, where the grid moves with the 
velocity of the stellar fluid, the latter reads: 

/ dF(DtS)=Dt/ dVE- I dV{EdWv). (17) 

J AV JaV JaV 

Here and in the following the symbol S represents one of 
the moments J or H (and also J' := er^J or H := e~^iJ). 



The computational domain 



< r < r„ 



is covered 



by Nr radial zones. As the zone center r^^i we define the 
volume- weighted mean of the interface coordinates and 
n+v 



1/3 



0, , 



1 . 



(18) 



Scalar quantities like ^i+i sue defined on zone cen- 

ters, whereas vectors like Hi and Li live on the interfaces 
of radial zones. "Miscentered" values like, e.g,. Ji or iJ^+i, 
are in general computed by linear interpolation, using the 
values of the nearest neighbouring coordinates (Ji_i and 
Ji+i, or Hi and i/^+i, respectively). 



3.3.2. Finite differencing of the moment equations 



Applying the procedures described in Sect. |3.3.1| to the 
moment equations (Eqs. 0, H) we obtain for each neutrino 
species ly G {vc-,v, 



. } (the corresponding index is 
suppressed in the notation below) the finite differenced 
moment equations. On an Eulerian grid they read: 

jn+l _ jn 



- — 



i+i 



[{I- Ik) 



dr 



f jn+l 



9/3 



dt 



[//f j"+Vi,+i=c(°):'^,,+i, (21) 



Spectral discretization; The spectral range < e < emax is 
covered by a discrete energy grid consisting of energy 
"bins" , where the centers of these zones are given in terms 
of the interface values as 



1 , 



1 



(19) 



and the width of an energy bin is denoted by Ae^^j. := 
Cj+i — tj - The radiation moments are interpreted as aver- 
age values within energy bins: 



1 



deS,(e). 



(20) 



Time differencing: We employ a time-implicit finite differ- 
encing scheme for the moment equations (Eqs. ^ H) in 
combination with the equations which describe the ex- 
change of internal energy and electron fraction with the 
stellar medium (Eqs. ^ ^6|). Doing so we avoid the re- 
strictive Courant Friedrichs Lewy (CFL) condition which 
explicit numerical schemes have to obey for stability rea- 
sons: For neutrino transport in the optically thin regime 
the CFL condition would require cAt < Ar, where Ar is 
the size of the smallest zone and At the numerical time 
step. Even more importantly, the time-implicit discretiza- 
tion allows one to efficiently cope with the different time 
scales of the problem: The "stiff" source terms introduce 
a characteristic time scale proportional to the mean free 
path A of the neutrinos, cAt = A = Ar/Ar, which can be 
orders of magnitude smaller than the CFL time step in the 
opaque interior of a protoneutron star, where the optical 
depth At of individual grid cells becomes very large. 



rrn+l TTn 



.UkJ 



n+ll 



n+i 



A( 



dp 



dr 



H 



d(3 



dt 



[(1 + M J 



n+l] 



(22) 



Note that by virtue of introducing the "Eddington fac- 
tors" (see the following section for their actual computa- 
tion) Jh := H/J,fK := K/J and Jl ■= L/J, the system 
of moment equations has been closed. The "flux" in energy 
space across the boundaries between energy zones appear- 
ing in the fourth and fifth line of Eq. (^ ) and Eq. (|2^), re- 
spectively, is given in terms of the interface values Jj+i j, 
which are defined as the geometrical mean of j_i 
and Ji-^-i j+ii ^-iid the advection velocities 



'^+h3 



dp 



dt 



dp' 



dr 



[Ik 



(23) 



i+k 
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dr 



[h 



dp 



at 



(24) 



The nonlinear interpolation of J in energy space turned 
out to be crucial for obtaining a numerically stable algo- 
rithm. The Eddington factors fn, Jk and are interpo- 
lated linearly in energy space in order to avoid problems 
in case the odd moments {H or L) change their signs. 

Since the radiation moments are defined on a staggered 
mesh, the finite-difference equations are second-order ac- 
curate in space. First-order accuracy in time is achieved 
by employing fully implicit time-differences ("backward 
Euler"). Only the advection terms in the third and forth 
line of Eq. (|2^) and Eq. ( p^ ) , respectively, are treated dif- 
ferently (see Sect. 4.3.1 for a motivation): 



:= (1-C)S"+CS' 



(25) 



For stability reasons, the interpolation parameter ^ = 0.51 
was chosen slightly larger than 0.5, which makes the treat- 
ment of the advection terms "almost" second-order accu- 



rate in time (see, e.g., Dorfi 1998). The interface values 
and are taken as the corresponding upwind quantities 
with 



i + 



for A > , 
else . 



(26) 



In order to damp spurious oscillations behind sharp 
neutrino pulses travelling through nearly transparent 
zones, an additional diffusive term is added to the left- 
hand side of the firs t moment equation (E q. ^2|) . Its finite- 
differenced form is (Blinnikov et al. li 



^visc — 



2r? 



H 



n+l -A _ Tjn^ 



n+l ^2 



n+l - Ti 



Tjn+l 



1 



(27) 



In effect this term resembles the diffusivity of the donor 
cell scheme and hence allows one to selectively switch 
from the second order radial discretization to a first order 



scheme. Following Blinnikov et al. (1998) the (energy de- 



pendent) "artificial radiative viscosity" A is set to unity 
in transparent zones and goes to zero for optical depth 
Ar > 1. By construction the diffusive term vanishes when 
the luminosity Hr^ does not depend on radius and hence 
is only active during transient phases when sharp neutino 
pulses propagate across the grid. 

Using a similar diffusive term in the zeroth moment 
equation (Eq. |2^) allowed us also to overcome stability 
problems with the numerical handling of the velocity- 
dependent terms in the general form of the radiation mo- 
mentum equation, Eq. (Bh. 



3.3.3. Boundary conditions 

For the solution of the moment equations (Eqs. 0, ||), 
boundary conditions have to be specified at r = rmin, 
T — fmax, e = and e — Cmax- In the radial domain the 
values of quantities at the boundaries are iteratively ob- 
tained from the solution of the Boltzmann equation (see 
Sect. 3.4.3 for the boundary conditions employed there), 
which has the advantage that in the moment equations no 
assumptions have to be made about the angular distribu- 
tion of the specific intensity at the boundaries. Specifically, 
at the inner boundary TJg i is given by Hit., r = r^-m, e) 
as computed from the Boltzmann equation. To handle the 
outer boundary, we ap ply Eq. (p^) on the " half shell" be- 
tween rj^^_i and r^v^ (Mihalas & Mihalas 1984). Where 
necessary, Jjy^+i j+i is replaced by /[/ff]jv^ j+i 

in Eqs. (|2l| |2|). 

At e = the flux in energy space vanishes and hence 
we simply set (oUi^i^QJ^_^l = in Eq. (^. At the up- 
per boundary of the spectrum the flux in energy space, 
is computed by a (geometrical) extrap- 



n+l 



J 

olation of the moments J,. 



and J,, 



and 



by a linear extrapolation of the advection velocity using 



for Eq. 



and ' 



H+i-N. 



_ 1 . Analogous expressions are used 



3.3.4. Finite differencing of the source term 

The finite differenced versions of the operator-splitted 
source terms in the energy and lepton number equations 
(Eqs. |l^, |l^) of the stellar fluid read: 

P^+l _ n 



fn+l _ fr, 



An 



n+l 

2 'J^ 2 > 



(28) 



47rm 



Pz+ 



2 J=0 



(29) 



3.3.5. Neutrino number transport 



In the flnite-differenced version of the (monochromatic) 
neutrino energy equation (Eq. pi] ) the derivative with re- 
spect to energy, d/de, has been written in a conservative 
form. When a summation over all energy bins is performed 
in Eq. ( ^l| ) , the terms containing fluxes across the bound- 
aries of the energy zones telescope and an appropriate 
finite differenced version of the total (i.e. spectrally inte- 
grated) neutrino energy equation is recovered. This essen- 
tial property does, however, not hold automatically also 
for the neutrino number density Af := Att/c dej^{e), 



when the naive definition J71 



i+i,]+hi' 



J, 



i+hd+h 



3 + h 
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is adopted. By inserting the latter expression into Eq. 
and summing over all energy bins, it can easily be verified 
that the corresponding fluxes across the boundaries of the 
energy bins do not cancel anymore due to the finite cell 
size Ae,- 1 i . 

In order to avoid this problem, we instead derive the 
moment equations for J and Ti, (by multiplying Eqs. ^ || 
with €~^) and recast them into a conservative form: 



,ld_ 

'cdt 



dr 
d_ 
d'e 



dr 

r or c ot 



dr 



dt 



(30) 



,ld_ 

'cdt 



d 



dr 
d_ 
'd'e 



1 d 

[n-c) 



JC-J 



dp ^ Idp 



di 



-C 



dt 



K 



r dr 



c dt 



(31) 



Equations ( |30| , |3l| ) are then discretized as a separate set 
of equations. Consequently, the finite-difference represen- 
tation of the monochromatic number density J^^^i j^i 
and the monochromatic number flux 7i, , , i are now con- 
sidered as new variables, which are determined by solv- 
ing an extra set of two moment equations. Note that 
treating and Ti^ jj^i as additional variables is 

not necessarily in conflict with the relations J{e) = 
J(e)/e and 7^(e) = H{e)/e {JC and C are defined hke- 
wise), which constitute dependences between the quan- 
tities. Since J^+i j+i and J^+i j+i, and H.^+i 
are defined as average values within the energy interval 
[ej, Ej+i] of different energy moments of the same function 
(rather than values of a function at point ej+ 1 ; cf. Eq. pO[ ) 

they can be considered as variables that are indepen- 3.4.1. Model equations 



Newton- Raphson iteration (e.g Press et al, 1992), using 
the aforementioned boundary conditions. 

This requires the inversion of a block-pentadiagonal 
system with 2Nr + 1 rows of blocks. The dimension of 
an individual block of the Jacobian is (4iVe -I- 2) in case 
of the transport of electron neutrinos and antineutrinos 
(the number of energy bins is multiplied by a factor 
of two because and are treated as separate parti- 
cles, the other factor of two comes from the additional 
transport of neutrino number). In contrast, muon and tau 
neutrinos and their antiparticles are considered here as one 
neutrino species because their interactions with supernova 
matter are roughly the same. In the absence of the cor- 
responding charged leptons, muon and tau neutrinos are 
only produced (or annihilated) as uv pairs and hence do 
not transport lepton number. Therefore the blocks have a 
dimension of (A^e -I- 1) for these flavours. 

For setting up the Jacobian all partial derivatives with 
respect to the radiation moments can be calculated ana- 
lytically, whereas the derivatives with respect to electron 
fraction and specific internal energy are approximated by 
finite differences. 

3.4. Formal solution of the Boltzmann equation 

In order to provide the closure relations (the "variable 
Eddington factors") for the truncated system of moment 
equations, we have to solve the Boltzmann equation. Since 
the emissivity rj and the opacity x depend in general also 
on angular moments of X, this is actually a nonlinear, 
integro-differential problem. However, 77 and x become 
known functions of only the coordinates i, r, e and /i, 
when J and H are used as given solutions of the moment 
equations. This allows one to calculate a formal solution 
of the Boltzmann equation. 



dent of each other, provided that the inequalities ej < 
< ej+i with (e)i+ij+i = J,+ i,j-+i /J^+i j+i 
or (e)j — j^i I are fulfilled (we thank our 

referee, A. Mezzacappa, for stressing the importance of 
this point). In practice, we find that this is always well 
fulfilled at high optical depths and low velocities. The con- 
straint is violated only in the steeply falling, high-energy 
tail of the neutrino spectrum. Systematically inspecting 
results from simulations we find that a part of the spec- 
trum is affected where J,- 1 i w^i Ae,- , 1 , ,1 is less than one 
percent of the maximum value, which means that the en- 
ergy integral over this spectral tail accounts for a negligible 
fraction of the total neutrino (energy) density. 



3.3.6. Numerical solution 

With the Eddington factors fn, fx, fh and flow 
field (3 being given, the nonlinear system of equations 
, |3C|, IsTI, |28l E9I) is solved for the variables 



Making the change of variabl e s (cf. Yorke 1980 ; Mihalas 
|fc Mihala£||l984| ; |Kornei) |l99^ ; [Baschek et al.||l997[ ) 



(r, /x) (s ^r,p r^Jl - y?), Ai > , (32) 

and introducing the symmetric and antisymmetric aver- 
ages of the specific intensity 



j{t,s,p) := i(X(^)+X(-Ai)), 

Eq. (^ can be split into two equations for j and h, 

ID 9 , 

-J + —h 



cDr 



ds 
d_ 
'd'e 



(33) 
(34) 

(35) 



(Eqs. 21 



2,/3 ^d(3 . „^ , 



dr 



{'^i-|-5J+3''^ij-|-7''^i-|-7j + 3''^i-|-5J>3''^i+5'-^ei+l} by & 



,1^ 

2d (3 



r dr c dt 
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cm 



ds' 
d_ 



(36) 



I dp 
c ot . 



+ 2 



9/3 
dr 



I dp 



The symmetric (se) and antisymmetric (ue) averages of 
the source term C are defined in analogy to Eqs. (|3^, Q). 
Following Burrows et al. (200C, with our Eq. ISSl correcting 



a misprint in their Eq. 24) we have made the replacements 



3^ / d(3 P\ dj 



- m') 



3 , dl3 p\ dh 



dr 



dfj, 



(3M^-l)(f --^-.(37) 
(V-2) (f -f)/.,(38) 



in deriving Eq. ( p5[) and Eq. (36), respectively, from 
Eq. (H). This modification can be motivated as follows: 
As a consequence of Eqs. (|3^, H^) no more partial deriva- 
tives of the intensity (and hence j and Ji) with respect to 
fj, appear in the model equations (Eqs. |3^, ^ ) of the origi- 
nal Boltzmann equation (Eq. This quality allows us to 
construct an efficient numerical solution scheme to be de- 
scribed in the next section. At the same time it is ensured 
that the monochromatic 0{f3) moment equations (Eqs. [7|, 
^) can be recovered when Eq. ( |35| ) and fi times Eq. ( ^q ) 
are integrated over angles. Thus, Eqs. (^, |3^) are consis- 
tent with the monochromatic moment equations without 



including the exact aberration terms (see also Burrows 
et al. |||OOOD . Mor cover the effects of angular aberration 
on the solution of the neutrino transport are found to be 
small (Liebendorfer, personal communication). Therefore 
we do not expect significant changes, in particular because 
only normalized moments H/J, K/J, L/J are computed 
from the solution of the Boltzmann equation to be used 
as closure relations in the moment equations. 

3.4.2. Finite-difference equations 

For the moment we shall ignore terms which contain fre- 
quency derivatives d/de (second line of Eqs. ^). These 
Doppler terms will be included in an operator-split man- 



the ray's intersections with the zone boundaries (centers) 
of the radial grid (cf. Fig. |^). With the "flux-like" vari- 
able h being defined at the zone boundaries Sik and the 
"density-like" variable j being defined at the zone cen- 
ters Sj^. 1 fc, the finite-differenced versions of Eqs. (35 , ^6|) 
finally can be written down (cf. Mihalas fc Klcinf 1982 , 
§V.2): 



J 



3 



A" 



I Ji+1 _ ,n+l 
_ „n+l 



ik+l,k 



,k + 5 



0, (39) 



ik,k 



ik,k 



J. 



ik+k.k 



^k'^^-k ik ^1^' 

+^::ljZ^k+Bt:kK::l-o^ m 

with the indices k = Kq, . . . , Nr — l and = fc, . . . , Nr — l. 
If an inner boundary condition at some finite radius r^i^ > 
is used, a number of "core rays" k = Kq, Kq + I, ■ ■ ■ , —1 
{—Kq G IN), which penetrate the innermost radial shell 
{pk < 7'min), Can bc defined in order to describe the angu- 
lar distribution of the radiation at the inner boundary. For 
details, see e.g., |Yorke| ( |l980| ), |KorneT| ( |1992| ), |Dorfi| ( |1998D . 
The coefficients A and B are combinations of the veloc- 
ity and angle-dependent terms and the source terms of 
Eqs. (|3^, ^). The definition of j and h will be given later, 
for the moment one may ignore the tilde symbol. We add 
a diffusive term to Eq. (Eol) analogous to the discussion in 
Sect. |3.3.2|, replacing Hr^in Eq. (|2^) by h. 



The frequency derivatives of Eqs. ( |35| , p6| ), which were 
ignored in the finite-difference versions, Eqs. ( p9| , ^ ), are 
taken into account in a separate step by operator-splitting 
(the partially updated values for j and h were marked by 
asterisks in Eqs. |4^). The discretization of the cor- 
responding terms is explicit in time and — provided the 
acceleration terms in the second line of Eqs. (35, ^6|) is 
neglected — can be performed straightforwardly using up- 
wind differences in energy space {I — 0, . . . , — 1): 



Equations (^, are then discretized on a so- 
called "tangent ray grid" (for an illustration, see Fig. |l|), 
the geometry of which being an immediate consequence 
of the transformation of variables given by Eq. (^^. 
Applying this transformation, partial derivatives of only 
one momentum-space coordinate s remain, whereas the 
second coordinate p appears only in a parametric way (cf. 
Baschek et al. 1997 ). This greatly facilitates the numerical 
solution of the system. 

A "tangent ray" k is defined by its "impact parame- 
ter" Pk — rk Sit ^ — 0. The coordinate s serves to measure 
the path length along the ray. On each tangent ray k, a 
staggered numerical mesh is introduced for the coordinate 
s. The zone boundaries (centers) of this mesh are given by 



with 



■ *n+l 



Ct 



n+1 



.<^i+iv..,.3 .,.,,(1+1) -^1^.,.] .,.,,(1) 



„4-n+l _ j-n 

i,n+l _ 7 *n+l I 

111 — ''■ ;^ 1 ^ 1 



I + ^ else , 



(41) 



J ; - i for V.,. > , 



(43) 
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and 



fl 







'd(3' 


r 


ik 


dr 



(44) 



where, for better readability, dashes replace the appropri- 
ate index pairs corresponding to radius and impact pa- 
rameter, {ik + 1/2, k) in Eq. (|l]) and {ik,k) in Eq. (||), 
respectively. 

It should be possible to proceed along similar lines for 
including the exact aberration effects in the Boltzmann so- 
lution which, for the time being, includes aberration only 
in an approximate way on the basis of the replacements 
given by Eqs. 



(37 



L In practice, however, the adopted 
tangent-ray geometry does not allow for a conservative 
discretization of the angular derivatives in a way as sim- 
ple as in the case of the frequency derivatives. We plan 
to spend extra work on the omitted angular derivatives of 
the intensity in order to include them in future version of 
our code. 

3.4.3. Boundary conditions 

On the tangent ray grid boundary conditions must be 
specified for each ray k at Sk^k and at SN^^k- At the inner 
core radius (i = 0; fc = Kq, . . . , and s_ i j, :— so.fc) we set 



with T, 



0,k 



the boundary condition j_ 
2^(^, ''minj Mo.fc)- For the remaining rays (fc — 1, . . . , N^), 
symmetry and Eq. (^) imply hk^k = 0, since iik,k — 0. 

At the outer radius {i = Nr\ k = Kq, . . . , Nr and 
s'Zi^ k ■= *A^,fc) consider Eq. (|^) on the "half shell" 
between r 



i and rjVr (cf. Sect 3.3.3) and make the re 



placement ijv^+^ 



+hNr k, withX 



The physical boundary conditions are described by the 
functions 2{t, Tmin, /i) and X(t, Tmax, — /^) with < < 1, 
which specify the specific intensity entering the computa- 
tional volume at the inner and outer surfaces, respectively. 

At the boundaries of the energy grid we use the same 
type of boundary conditions as described for the moment 
equations (cf. Sect. [3.3. 3 ). 



3.4.4. Numerical solution 

By virtue of the approximations used to derive the model 
equations (Eqs. |36| ) and upon introducing the tangent 
ray grid, the system of Eqs. (39, 4^) with suitably chosen 
boundary conditions can be solved independently for each 
impact parameter pk, each type of neutrino (note that for 
simplicity we have dropped the index u in our notation), 
and — because Doppler shift terms are split off — each 
neutrino energy bin e^^i (index also suppressed). 

For the same reasons as detailed in Sect. 3.3, we have 



employed fully implicit ( "backward Euler" ) time differenc- 
ing. Solving Eqs. (^, pO| ) therefore requires the separate 
solution of Nk x x {Nk Nr - Kq + 1 is the 
number of tangent rays) pentadiagonal linear systems of 
dimension < Nr. On vector computers, this can be done 



very efhciently by employing a vectorization over the index 
k. Once the numerical solution for j and h has been ob- 
tained, the monochromatic angular moments and thus the 
Eddington factors fn = H/J, fx = K/J and /l = L/J 
can be computed using the numerical quadrature formulae 



J{n) 

K{n)-- 
Un) = 



k=Ko 



pi i 

I dp.fih{n,p,) « ^ hkhk 

^0 k=Ko 



^ ^ jik^ik 
k=Ko 



- / dfif?h{ri,fi)x ^ hikd 

^0 k=Ko 



(45) 
(46) 
(47) 
(48) 



Explicit expressions for the quadrature weights Uik , bik 
and Cik can be found in Yorke (198C), dik is calculated 
in analogy. 

Unless the velocity field vanishes identically (in this 
case j = j and h = h in Eqs. ^ and |4^, respec- 
tively), an additional complication arises, since the par- 
tial derivative D/Dt has to be evaluated not only at fixed 
Lagrangian radius (the index i in our case) but also for 
a fixed value of the angle cosine fi. The angular grid 
{tA!^k}k=Ko...i, which is associated with the radius r", how- 
ever, changes as moves with time. As a consequence, 
e.g., the values of h at the old time level used in Eq. (40) 
are h^k := /i(t" 



^r/^) and not h. 



ik 



hit"" 



the solution of the previous time step (cf. Mihalas fc 
Mihalas 1984, §98). At the beginning of the time step 



we therefore have to interpolate the radiation field for 
each fixed radial index i from the {p-"k}k=Ko...i-S'-'^d onto 
the {/i"^+^}fe=if(,...i-grid. If an Eulerian grid is used for the 
computation, an additional interpolation in the radial di- 
rection is performed to map from the fixed r^-grid to the 
coordinates given by r"'"^ :— ri — ViAt (see Fig. |^). This 
seemingly complicated procedure has the advantage of be- 
ing computationally much more efficient than directly dis- 
cretizing Eqs. ( ^ , ^6| ) in Eulerian coordinates. In this case 
the vdr-teim, which originates from making the replace- 
ment Dt = dt + vdr , would couple grid points of different 
impact parameters and would therefore defy an indepen- 
dent treatment of the tangent rays. 



3.5. Iteration procedure 

A transport time step starts by adopting guess values (e.g., 
given by the solution of the previous time step) for J, H, 
Yc, and e. The quantities 1^0, e and the density p deter- 
mine the thermodynamic state and thus the temperature 
and chemical potentials. These, together with J and H 
are used to evaluate the rhs. of the Boltzmann equation. 
This allows one to compute its formal solution and the 



Eddington factors as detailed in Sect. 3.4. The Eddington 
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Fig. 1. Eulerian radial and tangent ray grid of the com- 
putation (solid lines) and virtual Lagrangian shell at time 
level = t^~^^ — At (dotted). For simplicity we have 
assumed a radial motion of only one shell with index i. 
The filled circles mark the positions where j and h are 
known as the solutions of the previous time step on the 
Eulerian grid. For v ^ these positions do not coincide 
with the locations (?'°''^, Mi/c"*^^): where initial values of the 
radiation field are required for computing the current time 
step (open circles). The latter positions are subject to the 
requirement that the D/Dt operator has to be evaluated 
at fixed radial Lagrangian and angular coordinates (The 
radial dashed lines connect points with constant angle 
cosines fj, of tangent rays.). 



factors are then fed into the system of moment and source- 
term equations, the solution of which (Sect. yields 
improved values for J, Yc, and e. At this point, where 
required (e.g., for evaluating C^"^), the neutrino number 
density and number flux density Ti are conventionally 
replaced by J and e^^H. This procedure is iterated un- 
til numerical convergence is achieved. With the Eddington 
factors being known from the described iteration proce- 
dure, the complete system of source-term and moment 
equations for neutrino energy and number is solved (once) 
in order to accomplish lepton number conservation. In 
this step a nd H are treated as additional variables 
(cf. Sect. |3.3.5|) . 



3.6. Details of coupling neutrino transport and 
hydrodynamics 

Our experience with the operator-splitting technique has 
shown that considerable care is necessary in how pre- 
cisely the equations of hydrodynamics are coupled with 
the neutrino transport and how the fractional time steps 
are scheduled. In the following we therefore describe the 
used procedures in some detail. 



Since the hydrodynamics code and the transport solver 
in general have different requirements for numerical res- 
olution, accuracy and stability, the discretization in both 
space and time of the two code components is preferably 
chosen to be independent of each other. In supernova sim- 
ulations, for example, the time step limit (from the CFL 
condition) for the hydrodynamic evolution computed by a 
time-explicit method is typically more restrictive than in 
the time-implicit treatment of the transport part. Since 
the transport is computationally much more expensive, 
one wants to use a larger time step than for the hydro- 
dynamics part. This means that a transport time step is 
divided into a suitable number of hydrodynamical sub- 
steps. 

3.6.1. Schedule of updates 

Starting from the hydrodynamic state (p", e", i^", d") at 
the old time level f\ the time-explicit PPM algorithm first 
computes a solution of the hydrodynamical conservation 
laws without the effects of gravity and neutrinos. From 
the updated density, the gravitational potential can be 
calculated by virtue of Poisson's equation and finally the 
gravitational source terms are applied to the total energy 
equation (Eq. ||) and the momentum equation (Eq. 
This completes the "hydrodynamics" time step with size 
A^: Hvd, which is limited by the CFL-condition (see e.g., 



LeVequq |1992D . By performing a total of N^^^ substeps 
the hydrodynamic state is evolved from the time level 
to the new time level given by 



Hyd • 



(49) 



i=l 



where N^^^ e IN is determined by the condition At" < 
At^^^. Here /S.t^^^ is the maximum size for the transport 
time step, which is compatible with the constraints that 
the relative changes are at most 10% for the monochro- 
matic neutrino energy and number density, and the rela- 
tive changes of the density, temperature, internal energy 
density and electron fraction must not exceed values of 
a few percent during the time interval At". This leads to 
a partially updated hydrodynamic state (p"^^, e*, Y^* , v*) 
which includes the effects due to hydrodynamic fluid mo- 
tions and the acceleration by gravitational forces and neu- 
trino pressure (see below). The corresponding state vari- 
ables are then mapped onto the transport grid. During the 
transport time step of size At", the transport equations 
in combination with the evolution equations for the elec- 
tron fraction and internal energy of the stellar medium 
in response to neutrino absorption and emission (Eqs. pS] , 
p9| ) are solved. This yields the source terms for energy 
and lepton number as Qe — p"^^ ■ (e"+^ — e*)/At" and 
= /3"+^ • (Fc""^^ - yc*)/At". Since the transport grid 
in general does not necessarily coincide with the hydro- 
dynamics grid, the new specific internal energy e"^^ and 
electron fraction Yc"^~^^ as computed on the transport grid 
do not directly describe the new hydrodynamic state. 
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Instead of interpolating the quantities e"'*'^ and Fo"^^ 
back onto the hydrodynamics grid, we map the source 
terms by a conservative procedure and then update the 
electron number density and the total energy density on 
the hydrodynamics grid according to 



=p"+ie* +(Q^+i + v'^+'Qlt') ■ Ai" , (50) 
[pY,r+' =p"+ii;*+QJ^+i . A<" . (51) 

Equations ( |50| ) and ( |5l| ) express the effective influence of 
the source terms Qe + vQm and Qn over the time of a 
transport step, but in fact we apply the corresponding 
sources (as given at the old time level) during each sub- 
step of the hydrodynamics. The quantities e* and Yc* as 
needed for the transport part of the code, are recovered 
a posteriori by subtracting the accumulated effects of the 
neutrino sources again. The momentum equation of the 
stellar fluid is treated in a similar way. Accounting for 
the momentum transfer (acceleration) by neutrinos after 
each individual hydrodynamical substep (by using the old 
momentum source term Q^), however, is of crucial im- 
portance here. When a sizeable contribution of the neu- 
trino pressure is ignored during a larger number of hydro- 
dynamical substeps a potentially hydrostatic configura- 
tion can be severely driven out of mechanical equilibrium. 
Correspondingly, the fluid velocity v* used in the trans- 
port step includes the effects due to the acceleration by 
neutrinos. 

After the transport time step has been completed, the 
new neutrino stress is used for correcting v* to give 

the new velocity v"'^^ at time level i""*"^: 



[pvr+' = p 



(Qm 



n+1 



- Qm) ■ At" . 



(52) 



3.6.2. Arrangement of spatial grids 



Radial discretization of the transport and hydrodynamic 
equations is done on independent grids. Therefore there is 
freedom to choose the number of zones and their coordi- 
nate values separately in both parts of the code. 

Only quantities which obey conservation laws have to 
be communicated from the hydro to the transport grid. By 
invoking the equation of state all other thermodynamical 
quantities can be derived from the density, the total energy 
density, momentum density and the number densities of 
electrons and nuclear species. In the other direction it is 
the neutrino source terms which have to be mapped back 
onto the hydro grid. 

For the mapping procedure we assume the conserved 
quantities to be piecewise linear functions of the radius in- 
side the grid cells, with parameters determined by the cell 
average value s and so-called "monotonized slopes" (for de- 
tails see, e.g., Ruffert 1992 , and references cited therein). 
In order to achieve global conservation of integral values 
we then simply average this function within each cell of the 
target grid of the mapping (cf. Dai & Woodward 1996). 



Using this procedure in dynamical supernova simula- 
tions, we discovered spurious radial and temporal fluctu- 
ations of the temperature, electron fraction, entropy and 



related quantities inside the opaque protoneutron star un- 
less the hydro and transport grids coincide exactly in this 
region. The scale of these fluctuations is sensitive to the 
radial cell size and the size of a time step, the relative 
amplitude was found to be typically of the order of a few 
percent, with the exact number varying between differ- 
ent quantities (see Rampp & Janka 2000| ). These fluctu- 
ations are understood from the fact that the mapping of 
the source terms on the one hand and the internal energy 
density (resp. temperature) on the other hand implies de- 
viations from local thermodynamical equilibrium between 
neutrinos and stellar medium. The attempt of the neutrino 
transport to restore this equilibrium within a given grid 
cell leads to large net production or absorption rates, driv- 
ing the temperature in the opposite direction and causing 
it to perform oscillations in time around a mean value. 
Despite these oscillations and fluctuations, the use of dif- 
ferent numerical grids inside the protoneutron star did not 
cause noticeable problems for the accuracy of the "global" 
evolution of our models, because the conservative mapping 
describes the exchange of energy and lepton number be- 
tween neutrinos and the stellar fluid correctly on average. 

Being cautious, we take advantage of the option of us- 
ing different hydro and transport grids only during the 
early phases of the collapse and in regions of the star, 
where neutrinos are far from reaching equilibrium with 
the stellar fluid. In this case no visible fluctuations occur. 



3.6.3. Conservation of total energy and lepton number 

The numerical treatment of the radiation hydrodynamics 
problem should guarantee that the conservation laws for 
energy and lepton number are fulflUed. 

Provided the acceleration terms (oc d(3/dt), which are 
of order {v^ / c^) and thus usually very small compared to 
the other terms, are ignored, the constancy of the total lep- 
ton number is ensured by, (a) a conservative discretization 
of the neutrino number equation (Eq. ^o|) , (b) a conser- 
vative handling of the electron number equation (Eq. 
and (c) the exact numerical balance of the source terms 
(cf. Eq. H) -47rmB/dy /g°°deC(°)(e) (defined on the 
transport grid) and J dV Qn (defined on the hydro grid) . 
Point (a) requires that in Eq. ( |30| ) the fiux divergence 
is discretized in analogy to the second line in Eq. (21) 
and that the jSdJ/dr and {2j3/r + dl3/dr)J terms are 
combined to d\v{(3 J) to be discretized in analogy to the 
third line in Eq. ( |2l| ) . The energy derivative in Eq. (^0|) i s 
treated in a conservative way as described in Sect. (3.3.5). 



Point (b) is achieved by the use of a conservative numeri- 
cal integration of the electron number equation (Eq. ^ in 
the spirit of the PROMETHEUS code, and requirement 
(c) is fulfilled by employing a conservative procedure for 
mapping the electron number source term from the trans- 



port grid to the hydro grid (see Sects. 3.6.1, 3.6.2). Doing 
so, the total lepton number remains constant in principle 
at the level of machine accuracy. 
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Different from the number transport, where the zeroth 
order moment equation for neutrinos by itself defines a 
conservation law, the derivation of a conservation law for 
the total energy implies a combination of the radiation 
energy and momentum equations. The use of a staggered 
radial mesh for discretizing the latter equations defies a 
suitable contraction of terms in analogy to the analytic 
case. Therefore our numerical description does not con- 
serve neutrino energy with the same accuracy as neutrino 
number and the quality of total energy conservation has to 
be verified empirically for a given problem and numerical 
resolution. 

For our supernova simulations, tests showed that neu- 
trino number is conserved to an accuracy of better than 
10^^^ per time step, while for neutrino energy a value 
below 10^^ is achieved. With a typical number of about 
50 000 transport time steps for a supernova simulation we 
thus find an empirical upper limit for the violation of en- 
ergy conservation of 0.5% of the neutrino energy. This 
translates to 0.05% of the internal energy of the collapsed 
stellar core, i.e. a few times 10'^^ erg in absolute number. 
Errors of the same magnitude are introduced by the non- 
conservative treatment of the gravitational potential as 
a source term in the fluid-energy equation (Eq. ^). Note 
that the use of different grids for the hydrodynamics and 
the transport does not affect the energy budget because 
we employ a conservative mapping of the neutrino source 
term between the grids (see Sects. 3.6.1, 3.6.2| ). 



and the neutrinos: 




PtotC 



(53) 



where ptotc^ p{c^ + e) is the total ( "relativistic" ) en- 
ergy density and P — Att/c de K the neutrino pres- 
sure. The calculation of the gravitational mass m(r) := 



/;dr'47rr"(ptot + c-^E 



^UF/T) takes into ac- 



count contributions of neutrino energy density E = 
Att/c /p de J and flux F — An deH. The metric func- 



tion r is calculated as T{r) = \/l + U{rY — 2Gm{r)/rc'^ 
with the term accounting for the effects of fluid motion. 

Equation (|53|) can be used in the Newtonian hydrody- 
namic equations (Eqs. |^, ^ in order to approximately take 
into account general relativistic effects (cf. Kcil 1997 ). The 
quality of this approach has to be ascertained empirically 
by comparison with fully general relativistic calculations. 
In our case such a comparison yields quite satisfactory 
results (see Sect. 4.3). 



3.7. Approximate general relativistic treatment 



3.7.2. Approximate GR transport 

The general relativistic moment equations describing 
transport of neutrino energy, momentum and neutrino 
number can b e der ived from the Lindquist-equation 
(cf. Eq. I, Sect. |2.2.l[ ). They are: 



We have not yet coupled our general relativistic version 
of the neutrino transport to a general relativistic hydro- 
dynamics code. For the time being we work with a ba- 
sically Newtonian code, which was extended to include 
post-Newtonian corrections of the gravitational potential. 
We hope that the deeper gravitational potential can ac- 
count for the main effects of general relativity on stellar 
core collapse and the formation of neutron stars which 
do not approach gravitational instability to become black 



holes (cf. Brucnn et al. 2001). Because the general rela- 
tivistic changes of the space-time metric are ignored, a 
consistent description of the neutrino transport requires 
that the fully relativistic equations are simplified such that 
only the effects of gravitational redshift and time dilation 
are retained. 



cDt 



i__^(i?2^C*) 



d_ 

We 



+ e'* 



.U 
R 



4(. 



{3J-K) 



ra^e* H 

c-^DtAK 



K) + c-^DtA K + TdRB^ 
-iDtA(J + i^) = e* C^°\ 



1 D 



r d 



de 



(4 



r^flc* J + e*r 



cr^BtAL 



R 



(54) 



-iDtA)i/ = e*C(i), (55) 



3.7.1. Modified gravitational potential 

By comparing the Tolman-Oppenheimer-Volkoff equa- 
tion for hydrostatic equilibrium in general relativity (see, 
e.g., Kippcnhahn & Wcigcrt 1990, Sect. 2.6) with its 
Newtonian counterpart (cf. Eq. ^ ane can define a mod- 
ified "gravitational potential" which includes correction 
terms due to pressure and energy of the stellar medium 



for the energy transport, and 



1 D 

cDt 



d_ 



r d 



i?2 dR 

u 

R 



[R^ne^) 



K.) 



TdBe^ n 



R 



(2e*-5-f c~iDtA)^ = e*C(*", (56) 
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d_ 
We 



V R, 



+ e'^^{n + C)+ c-^DtA {2n - £) = e* C^^^ , (57) 

for the number transport. 

In the approximate treatment we neglect the distinc- 
tion between coordinate radius and proper radius {Or -^■ 
dr, r = 1) in Eqs. ([54|-p7|), and identify corresponding 
quantities with their Newtonian counterparts (R — > r, 
e* [/ ^ P, c~^DtA — > d/3/dr). The same approximations 
are made in the "parent" Boltzmann equation from which 
the moment equations of the relativistic approximation 
can be consistently derived. 

Accordingly, the approximations to Eqs. (p^-|57|) con- 
tain only general relativistic redshift and time dilation 
effects for neutrinos. Coupling the transport with the 
Newtonian equations of hydrodynamics, these restrictions 
to a fully relativistic treatment are necessary in order to 
verify conservation laws for energy and lepton number of 
the coupled system. 

Finite-difference versions of the moment equations and 
the co rresponding parent Boltzmann equation for the ap- 



global deformation (e.g., due to rotation) in layers which 
are opaque to neutrinos, but where inhomogeneities and 
anisotropies are present only on smaller scales (e.g., due 
to convection). Multi-dimensional hydrodynamical simu- 
lations suggest that convective processes occur in two dis- 
tinct regions of the supernova core: 



(a) 



Convection inside the opaque protoneutron star causes 
deviations of the structural and thermodynamical 
quantities (like density or temperature) from spheri- 
cal symmetry typically of the order of a few percent 



(b) 



(Keil et al. 199(:; Keil 1997). Moreover, the time scale 
for changes of such local fluctuations is short compared 
to the neutrino diffusion time scale. Hence, on the neu- 
trino diffusion time scale no persistent local gradients 
in lateral and azimuthal directions are present in the 
dense stellar interior. This allows one to disregard the 
neutrino transport in these directions in a first approx- 
imation, which is at least correct in the temporal av- 
erage. 

Convective overturn motions occur between the neu- 
trinospheres and the stalled hydrodynamic shock. 
Deviations from spherical symmetry are present there 



Burrows et al 



techniques described in Sects. 3.3 and 3.4. The lapse tunc- 



proxinate GR transport are obtained by applying the Janka 1997 ; Mezzacappa et al. 1998b ; Kifonidis et al 



on larger angular scales (Herant et al. 1992, 1994 
1995|; iJanka & MiiUeil |1996|; Huller & 



tion is calculated by integrating the general relativistic 
Euler equation din /dr — —{ptotc^ + p)^^ {dp/dr — 
Qm/Y) inward from the s urface, where th e boundary con- 
dition e* = r is applied ([van Riper 1979| ). 



20001 ), but the neutrinos are only loosely coupled to 
the stellar medium in these regions. Therefore neutrino 
transport in the lateral directions does not play an im- 
portant role and cannot have a significant influence on 
the local emission and absorption of neutrinos. 



3.8. Approximate multi-dimensional neutrino transport 3.8.2. "Ray-by-ray" transport 



Multi-dimensional frequency dependent neutrino trans- 
port in moving media and relativistic environments is a 
challenging problem for future work. Since convective phe- 
nomen a were recognized to be highly important in super- 
novae dHerant et al.| |1994 furrows et al.| |l995| ; [janka fc| 
Miiller|fl99e| ; |Keil et al.||l996i and rcfs. therein) one wouk 
of course like to perform simulations with Boltzmann neu- 
trino transport also in two and three dimensions. Here 
we suggest an approximate approach based on a straight- 
forward generalization of our variable Eddington factor 
method, which offers some advantages concerning com- 
putational efficiency. The approximation should be con- 
sidered as an intermediate step between spherically sym- 
metric and fully multi-dimensional models. Of course, 
the quality of the approximation for the neutrino trans- 
port which we will describe below, will finally have to 
be checked by a comparison with fully multi-dimensional 
transport calculations. 

3.8.1. Basic considerations 

Our approximate treatment may be a reasonably accurate 
and physically justifyable approach for describing neutrino 
transport in situations where the star does not show a 



Under these circumstances the specific intensity 
T{t, r, "(9, if, e, n) can be assumed to depend mainly 
on r but only weakly on longitude tp and latitude 1} of 
the background medium. Hence, like in the spherically 
symmetric case, the dependence of the specific intensity 
on the direction of propagation n can be described 
by only one angle fj, := n ■ r/\r\. The flux is thus 
approximated as = 47rr/|r| • H{t,r,'d,ip,e) and the 
scalar P = Air/c ■ K{t, r, (p, e) is sufficient to define the 
radiation stress tensor. 

In the moment equations, gradients in the and the 
(/9-direction, which describe the transport of energy and 
neutrino number into these lateral and azimuthal direc- 
tions, are neglected, yet the parametric dependence of the 
(scalar) moments on the coordinates i9 and (p is retained 
and the radial transport is computed using the moment 
equations independently in each angular zone of the stel- 
lar model. In order to close the set of moment equations, 
variable Eddington factors are computed by iterating the 
Boltzmann equation and the corresponding moment equa- 
tions on a spherically symmetric "image" of the stellar 
background. The latter is defined as the angular average 
of physical quantities ^ € {p,T,Ye, (3, . . .} according to 
C(t,r) := (47r)-i/J^i^dcosi? J^'' dip ^{t,r,§, ^p). Note that 
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the variable Eddington factors are normahzed moments of 
the neutrino phase space distribution function and thus 
should not show significant variation with the angular 
coordinates of the star. This justifies replacing them by 
quantities that depend only on the radial position and 
time. 

Since for each latitude and longitude the moment 
equations (Eqs. 0, ^ in our approach are solved together 
with the evolution equations of electron fraction and in- 
ternal energy due to neutrino sources (Eqs. |l|, |l|), local 
radiative equilibrium can be attained properly and conser- 
vation of energy can be fulfilled. It is not obvious to us how 
these fundamental requirements could be met in a more 
simple approximation where one uses a one-dimensional 
transport scheme to compute the transport on a spheri- 
cally symmetric "mean star" , which is obtained at each 
time step by averaging the multi-dimensional hydrody- 
namical stellar model over angles. 

3.8.3. Computational efficiency 

Besides having significant advantages for easy and ef- 
ficient implementation on vector and parallel computer 
architectures, the suggested approach also saves com- 
puter time compared to a multiple application of a one- 
dimensional Boltzmann solver. Let Atfj^^ be the compu- 
tation time required for calculating a time step of the full 
one-dimensional transport problem and let n^* € IN be 
the number of steps that are necessary for the iteration 
between the moment equations and the Boltzmann equa- 
tion to achieve convergence. If n^_^ := • is the total 
number of angular zones that discretize the background 
star, the computation time for calculating one time step 
of the multi-dimensional problem with the described ap- 
proximate method is 



CPU 



1 



V \ A^CPU 



(58) 



Note that the iteration procedure (n^* > 1) has to be per- 
formed only once, namely for calculating the Eddington 
factors on the angularly averaged "image" of the back- 
ground medium. With typical values of n^* = 3 ... 10 and 
large values of n^,if, (~ 100) the computational effort for 
solving the multi-dimensional problem becomes about 



CPU 



(0.1. 



. 0.3) n^,^ • Aifi^u 



(59) 



n^.ip ■ Affp^, which would 



to be compared with Ai^^^ 
be necessary for a straightforward multiple application of 
the one-dimensional method to a number of n§^^ angular 
zones. 

4. Test cases 

Stationary solutions of the Boltzmann equation of ra- 
diative transfer can be calculated analytically, scmi- 
analytically or numerically to high accuracy in a few cases, 
where a particularly simple form of the emissivity and the 



opacity is assumed. These problems provide test cases for 
our new code for solving the equations of neutrino trans- 
port. To this end we impose suitably chosen initial and 
boundary conditions and solve the time-dependent neu- 
trino transport equations numerically, until a stationary 
state is reached. The numerical solutions can then be com- 
pared to the reference solutions. This first class of prob- 
lems, presented in Sects. 4.1 and O, does neither refer to 



specific properties of neutrinos nor will the equations of 
neutrino radiation hydrodynamics, i.e. the coupling of the 
neutrino transport equation to the equations of hydrody- 
namics, be tested. Time-dependent hydrodynamical tests 
will be described in Sect. 



4.3 



4.1. Radiative transfer on a static background 

We first consider the simplest class of radiative transfer 
problems, namely those where the "background medium" 
is assumed to be time-independent and static, i.e. the ra- 
dial profiles of the emissivity and opacity do not change 
with time and the velocity and acceleration vanish every- 
where and at all times. As usual, the background medium 
is assumed to be spherically symmetric. 



Propagation of a light pulse: Here, our method for com- 
puting the variable Eddington factor (see Sect. ^^) by 
solving the Boltzmann equation with a finite difference 



scheme ( Mihalas fc Klein 1982 ) is compared with results 



obtained from a general quadrature solution along charac- 
teristics ( |Yorkc| |198(]| ; |K6rnei] |l992| ; [Baschek ct "aL| |l997D 
For convenience, we set c = 1 in this section. 

As a test problem we consider a central point source 
which emits radiation with an intensity Xq into a spher- 
ical, static and homogeneous shell of matter bounded by 
two spheres of radius tq and R > tq. The only interac- 
tion of the radiation with the medium is by absorption 
(C = — Kal). The central source is active during the time 
interval — tq <t< oo. Initially (i = 0), there is no radia- 
tion inside the shell. 

By symmetry, the intensity vanishes for all but the ra- 
dial direction of propagation. In the absence of scattering, 
it is therefore sufficient to follow the propagation of the 
light pulse along a single radial (characteristic) ray. The 
dependence of the intensity on the angle can thus conve- 
niently be suppressed in the notation, writing 



X+(t, s) := 5(/i- 1) •X(t,s,/i), 



(60) 



with s := r — ro and 6{x) being the Dirac J-function. The 
equations to be solved are (cf. Eq. ||) 



(61) 
(62) 



subject to the boundary conditions X+(f,0) = Xo (for 
t > Q) and T~ (t, Smax) — 0, and the initial condition 
X='=(0, s) = (for < s < Smax ■— R— To). The analytical 
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Fig. 2. Outwardly directed intensity X"*" (normalized to the true postfront value Tq) as a function of the position s, 
calculated with our Boltzmann solver at time t = 100 for different values of the computational time step. The bold 
line represents the analytical solution. Panel a shows Model "VAC" with no absorption, in Panel b Model "ABS" with 
the absorptive shell is displayed. Employing an additional diffusive term in the transport equations helps damping the 
spurious post-front oscillations, which are most prominent for At/As = 0.02 (cf. the solid line compared to the dotted 
line in the plot inserted in Panel a). 



(weak) solution can easily be verified to read: 
Xq • exp(— Ka s) for s <t, 



I+{t, s) 



and 







else. 



X-(i,s) = 0. 



(63) 



With our imp lementation of the finit e-difference method 
adopted from Mihalas fc Klein ( 1982 ), we computed two 



models with a shell spanning the range < s < 200. We 
use an equidistant radial grid with a resolution of As — 1. 
Model "VAC" assumes a vacuum shell (^a = 0), whereas 
the absorbing shell of model "ABS" is characterized by 
Ka = 0.01, resulting in a total optical depth of r = 2. 
Figure shows our results for X+ at the time t — 100, 



together with the analytical solution (Eq. 63). We define 
the position of the numerically broadened light front as 
the radius, where the average of the true pre and post- 
front value of the intensity is reached. In all cases the 
mean propagation speed is well reproduced, with some mi- 
nor loss of accuracy for very large time steps. The shape 
of the light front, however, deviates from the true solu- 
tion. One observes an artificial precursor together with a 
reduced intensity behind the front, both effects resulting 
in a spatial broadening of the front. A clear trend towards 
larger diffusive broadening with increasing time steps can 
be seen from Fig. |^ and Table ^ This phenomenon was 



also observed by Mihalas & Klein (1982), who reported 
very similar results for their tests. 



At/ As 


# 


model VAC 
5s (FD) Ss (CH) 


model ABS 
Ss (FD) Ss (CH) 


0.02 


5000 


6 


25 


6 


26 


0.1 


1000 


8 




8 




0.5 


200 


18 


16 


16 


20 


1.0 


100 


26 


1 


22 


1 


2.0 


50 


37 


26 


29 


26 


10.0 


10 


80 




55 





Table 1. Width of the light front Ss (defined as the num- 
ber of grid points for which X+ lies between 0.1 and 0.9 
times its true postfront value) as a function of the time 
step At (or the number of time steps # required for the 
light front to reach s — 100). Computations with th e fi- 
nite difference method ("FD" ; [Mihalas fc Klcin||l982D and 
with the method of characteristics ("CH"; Yorkc 1980) 
are compared. Numbers obtained with the latter method 
are taken directly from Yorke (198C, Tables 1,2). They 



were, in addition, confirmed by an implementation of this 



numerical scheme by Korner (1992) 



Table I summarizes our results for the width of the 
front and compares them to simulations done with an 



implementation of the method of characteristics (Yorke 
1980; Korner 1992). The latter method has the advan- 
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tage to exactly reproduce the analytical solution without 
any smearing of the front, if the light front is advanced 
by exactly one grid cell per time step. Though the abil- 
ity to reproduce the exact solution appears to be a very 
appealing property of the method of characteristics, the 
necessary condition At = As can hardly ever be accom- 
plished in realistic simulations. This is simply due to the 
fact that in typical astrophysical simulations the radial 
resolution can be varying over several orders of magni- 
tude, whereas in most applications the global value of the 
time step is determined in some region of the grid. When 
doing radiation hydrodynamical simulations, for example, 
one is usually interested in resolving time scales different 
from those given by the speed of propagating light fronts. 
Accordingly, the time steps can be very different from the 
optimum for describing the light front propagation. Our 
finite difference scheme therefore seems to be preferable 
to the method of characteristics as far as the resolution of 
travelling discontinuities is concerned (cf . Table |l|) . 

For time steps At significantly smaller than the ra- 



dial zo ne width, a case not discussed by [Mihalas fc Klein 
(|1982|) , we observe some spurious postfront oscillations 



(see also [Mihalas fc Weaver] |1982| , §5). Their maximum 



amplitude, however, does not grow with time, which is 
in agreement with a local von Neumann stability analy- 
sis for the fully implicit implementation of the equations 
(Mihalas fc Klein 1982). If present in practical applica- 



tions, the oscillations can be damped by employing an 
additional diffusive term in Eq. (|6^), which is active in 
nearly transparent regions (r < 1) of the computational 
grid (cf. Sects. |3.3.2| , |3.4.2| ; [Blinnikov ct al.| |1998D . Due 
to the correspondingly increased diffusivity of the finite- 
difference scheme the light front is smeared over a larger 
number of zones {Ss = 18 for At/ As = 0.02; see also 
the inserted pa nel in Fig, p^) as compared to the orig- 
inal method of Mihalas & Klein ( [1982 ). Note, however, 
that the representation of the front is still better than the 
results obtained with the method of characteristics. 



Homogeneous sphere: The so-called "homogeneous 
sphere" is a test problem frequently employed for radia- 



tive transfer calculations (e. 




Bruenn 


1985 




Schinder 


fc Bludman 


1989; 3mit et al. 


1997 


). Physically, 


one can 



think of a static, homogeneous and isothermal sphere of 
radius R radiating into vacuum. Inside the sphere, the 
only interactions of the radiation with the background 
medium are isotropic absorption and thermal emission of 
radiation. 

Despite of its simplification, the model has some im- 
portant numerical and physical properties which are typ- 
ically found in practical applications: Finite difference 
methods are challenged by the discontinuity at the sur- 
face of the sphere. Moreover, the sudden transition from 
radiation diffusion inside an optically thick sphere to free 
streaming in the ambient vacuum — a similar but less ex- 
treme situation arises in the neutrino-heating region of a 
core-collapse supernova — is a major source of inaccura- 



cies for approximate radiative transfer methods like, e.g., 
flux-limited diffusion. 

The model is defined by setting C — Ki,{b — T) with 
b = const., Ka = const, for < r < i?, and Ka = 6 = 
for r > R. Since for this choice of parameters, the right 
hand side of the Boltzmann equation does not contain 
terms that depend explicitly on the angular moments of 
the radiation field, the solution of Eq. (||), with /3 = 0, 
can be obtained analytically by computing a formal so- 
lution. With the boundary conditions T{r = 0,/i = 1) = 
I{r = 0,M = -1) and I{r = i?nuax,-l < M < 0) = 0, 
which account for isotropy (due to spherical symmetry) of 
the radiation field at the center and no incoming radia- 
tion at the outer boundary, respectively, the result for the 
stationary state is (see e.g., ^mit et aL 1997): 



X(r,M) 
Here s is defined as 



6(1 



(64) 



r^i + Rg{r,fi) if r < i?, -1 < ^ < +1, 
s(r,/x) := <( 2i?g(r,Ai) if r > i?, x < fj, < +1, 

else. 



with g(r,^l)■.= \|^-[-j (i-/i2) 



(65) 



and X := \ I 



R 



The solution given by Eqs. (|6J, 4T) depends on only three 
independent parameters, namely the radial position rel- 
ative to the radius of the sphere, r/R, its total optical 
depth T — KaR, and the equilibrium intensity b. The lat- 
ter merely acts as a scale factor of the solution. 

We employ a radial grid consisting of 213 radial zones 
to cover the range between r = and r = i?max ~ 12i?. 
Approximately 200 zones were distributed logarithmically 
between r « 0.0006i?max and r — i?max, about two thirds 
of which were spent to resolve the sphere. Initially we set 
X « everywhere and evolved the radiation field until a 
stationary state was reached. 

In Fig. ^ we display the stationary-state solutions for 
two models, one with t = b ~ A (Figs. ^,c) and another 
with T = 6 = 26 (Figs. ||b,d), representative of spheres 
with low and high optical depth, respectively. In general, 
our numerical results agree very well with the analytical 
solutions. The Eddington factor K/J (as calculated from 
the Boltzmann equation) and the flux factor H/ J (as ob- 
tained from the moment equations) both follow the exact 
values very closely (Figs. ||c,d), in case of the flux fac- 
tor over many orders of magnitude. In both models, the 
transition to free streaming is reproduced very accurately 
(Figs. ^,d). This region usually causes serious problems 



for flux-limited diffusion methods (see e.g., Bruenn 1985 
pmit et ah] 1997). Also the forward peaking of the radi- 



ation distribution at large distances from the surface of 
the sphere (cf. Eqs. p4, 4.1) is described excellently by 
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Fig. 3. Numerical (dotted and dash-dotted lines) and analytical (solid lines) stationary-state solutions vs. radius 
(normalized to the radius i?max of the outer boundary) of homogeneous sphere problems with "low" (r — 4, left 
column, Panels a, c) and "high" opacity (t = 26, right column. Panels b, d). Panels a and b show the zeroth moment 
of the intensity. In Panels c and d the Eddington factor Kj J (dash-dotted) and the flux factor 77/ J (dotted) are 
displayed. The inserted panels show the flux factor on a logarithmic scale. Thin vertical lines mark the surface R of 
the homogeneous sphere. 



our method: At r = i?max, the tangent ray grid yields 
approximately 150 angular grid points to resolve the radi- 
ation field from the central source, which has an opening 
half angle of only Q k, h° (cosO ~ 0.996). These num- 
bers should be compared with calculations employing the 
discrete angles (Sn) method: In time dependent neutrino 
transport calculations, typically less than 10 angular grid 



points can be afforded to equidistantly cover the range 
< cose < 1. 

Although we have used a geometrical radial zoning for 
our tangent-ray grid, we do not find any systematic ef- 
fects caused by a "bias" in angular binning as described 



by Burrows et al. (200C). Far from a central source, they 
report Eddington factors and flux factors that asymp- 
tote to between 0.96 and 0.98 instead of 1.0 in this case. 
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For example, we find values of K/J ~ 0.996602 for the 
Eddington factor and H / J = 0.999361 for the flux factor 
at r = i?max (in the model with r = 4), to be compared 
with K/J = 0.996636 and H/J = 0.998626 obtained from 
the analytical solution (Eqs. pil 4.1). 



In case of the mean intensity J we observe a system- 
atic trend towards larger deviations from the true solution 
for spheres with larger total optical depth r. In our "low 
opacity" case (r — 4, see Fig. ||a), there is hardly any dif- 
ference visible, whereas in the "high opacity" case (r = 26, 
see Fig. ^d) the numerical solution slightly overestimates 
the true solution in a narrow region beneath the surface of 
the sphere and underestimates it in the ambient vacuum 
region. A detailed analysis of the data yields values of 6% 
and 10% for the relative deviations from the analytical so- 
lutions in the "low opacity" and the "high opacity" case, 
respectively. 

This finding is caused by the steep opacity gradient 
near the surface of the radiating sphere. Since we keep 
our radial grid to be the same in all models, the transition 
from high optical depth (T(r) — {R — r) > 1, ii r < R) 
to transparency (T(r) = 0, if r > R) occurs in a bound- 
ary layer beneath the surface of the sphere, which is less 
well resolved when the opacity is very large. In the "high 
opacity" case, the semi-transparent layer is geometrically 
much thinner than in the "low opacity" case and therefore 
the radial gradients of the radiation density are steeper 
and would require more radial zones for a proper descrip- 
tion. In fact, in the "low opacity" case the radial zones of 
our grid near the surface of the sphere are optically thin 
{At « 0.15), whereas in the "high opacity" case the out- 
ermost zone has already an optical depth of At « 2.3. A 
better quality of the numerical results could be achieved 
by increasing the spatial resolution in the vicinity of large 
opacity or emissivity gradients and/or by using higher- 
order difference schemes. 

Due to the absence of scattering in the discussion 
of the homogeneous sphere problem, the solution of the 
Boltzmann equation does not require any information 
from the moment equations. No iteration between both 
parts of the code is necessary. Therefore this problem al- 
lows one to test the algorithm which solves the Boltzmann 
equation for the radiation intensity independently from 
the numerical solution of the moment equations. For the 
stationary state, we find that the radiation moments cal- 
culated by a numerical integration of the intensities are 
consistent with the moments directly obtained from the 
moment equations to within an accuracy of less than a 
percent. 



Static scattering atmospheres: Hummer & Rybicki (1971) 
published stationary-state solutions for the spherical ana- 
logue of the classical Milne problem. The model com- 
prises a static, spherically symmetric, pure scattering at- 
mosphere of some radius i?max, with a central point source 
that is emitting radiation isotropically with a constant lu- 
minosity Lq. Due to the presence of scattering, the prob- 



lem is of integro-differential nature and defies solution by 
simply computing the formal solution of the Boltzmann 
equation. 

The opacity of the atmosphere is assumed to be solely 
caused by isotropic scattering with a simple power-law 
dependence on radius: 



X{t,r) = Ks{r) = r 



<r < R„ 



(71 > 1) . (66) 



By a redefinition of the unit of length we have set a scale 
factor a > to unity. This factor appears in the orig- 
i nal fo rm Ks{r,t) — ar^"' used by Hummer fc Rybicki 



( |1971| , Eq. 1.1.). The emissivity vanishes within the at- 
mosphere (77 = 0) and no radiation is entering from outside 
(X(t, i?„iax, /i) = 0, for -l<n< 0). 

For a number of atmospheres defined by various com- 
binations of the parameters n g {1.5,2,3}, -R max G 
{0.1, . . . , 100} and Lq = {Anf, [Hummer fc Rybicki) (|l97l|) 



computed numerically the zeroth moment J and the 
Eddington factor K/ J of the stationary-state solution as 
a function of radius. They found values of better than one 
percent for the accuracy of their results. The stationary- 
state solution for the first moment H as a function of 
radius can easily be derived analytically from the zeroth 
order moment equation (Eq. 0, with /3 = 0): 



Hir) = 



La 



1 



(47r)2 



(67) 



This means that the luminosity is conserved {L(r) :— 
Anr^ ■ 47riJ(r) = ig)- [Hummer fc Rybickij ( [l97l| ) also gave 
useful asymptotic expressions applicable to regions of very 
low and very high optical depth, respectively. 

The idealized concept of a central point source with 
a given luminosity Lq is in practice modeled by impos- 
ing a suitable inner boundary condition at some finite 
radius i?rnin, which bounds the atmosphere from below. 
Since all atmospheres with a scattering opacity accord- 
ing to Eq. (^) become optically thick at sufficiently small 
radius, it is reasonable to employ the diffusion ansatz 



^{t, Rmin, m) — 2^0 + 



(68) 



for the radiation field at the inner boundary. In our pro- 
gram, only the quantity h{t, i?min, ^J■) = needs to be 
specified. Using Eq. (||) and H = fih (cf. Eq. H), 

the parameter Xi is easily verified to be Xi = 3i7 (i?niin) ~ 
3Lo/(47ri?,„i„)2. 

We evolved the transport to stationary-state solutions 
for the two sets of parameter combinations (n = 1.5, 
i?„nn = 0.001, i?„iax € {0.1,1,10,100}), and {n = 2, 
i?min = 0.01, i?max G {0.1,1,10,20}). The total optical 
depth at r = iJmin is larger than 55 for all models of 
the former class {n — 1.5) and larger than 90 for all of the 
models with n = 2. Hence, we verify that the inner bound- 
ary is placed at a radius which is sufficiently small to jus- 
tify the use of the diffusion approximation (Eq. |6^) for the 
inner boundary condition. In order to test the sensitivity 
of the results to the corresponding angular dependence of 
the intensity (Eq. p8|), we have calculated a model with a 
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Fig. 4. Stationary state solutions for selected radiation quantities as a function of radius obtained with our radiative 
transfer code (solid lines). Our results are compared with the reference solutions (crosses) and asymptotic solutions 
(dashed lines) of Hummer & Rybicki (1971). Panel a: mean intensity J (times r^) for the combination of model 
parameters {n = 2,Rjnin — 0.01,i?niax = 0.1); Panel b: mean intensity J for (n = 2,i?niin = 0.01,i?niax = 10); 
Panel c: r^J for models with (n = 3/2, i?min = 0.001, i?max € {0.1,1,10,100}). The vertical arrows indicate the 
radial positions, where r = 1 is reached in the different atmospheres; Panel d: Eddington factor K/J for the test 
case with a purely radial distribution of the radiation intensity at the inner boundary and the model parameters 
(n — 2,i?niin = 0.01,i?niax — 10). Note that the solid line connects values which are evaluated at the centers of the 
radial zones of our computational grid, the first one being located at ri — 0.05. 



different angular dependence: I(i, i?niin, /i) = (5(/^ — 1) -Xin. ated from below by the central point source. Figure ^ 

The parameter T^^ is chosen such that L(i?rriin) ~ Lq. This shows the Eddington factor of the stationary state for this 

prescription is appropriate for an atmosphere around a model. The numerical solution deviates visibly from the 

central hollow sphere with radius i?min, which is irradi- reference solution only in the innermost radial zones. 
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In the models with n = 1.5 the standard radial grid 
consisted of 200 logarithmically distributed zones with 
some additional zones giving higher resolution near the 
surface of the atmosphere, whereas 250 equidistant radial 
zones were used for the models with n — 2. Additional 
tests with 500 radial grid points were performed for 
some selected models without finding significant changes. 
Figure ^ shows some results of our tests together with 
the data of Hummer & Rybicki (1971). The agreement 



is excellent. In all models we were able to reproduce the 
analytical solution for the flux (Eq. |6^) with a relative 
accuracy of at least 10^"* throughout the atmosphere. 

[Kadiative] trans fer in general relativity: Schinder & 
Bludman ( 1989 ) considered the effects of a strong gravita- 
tional field in some of the radiative transfer problems dis- 
cussed above, in particular in the cases of a homogeneous 
sphere and a static scattering atmosphere. They computed 
stationary-state solutions numerically for the general rel- 
ativistic equations of radiative transfer and compared the 
results to the weak field limit. For a static background 
(9(i? = dtA = 0), and by applicatio n of a change of coor 
dinates (-R, e) i— > (i?, Coo ■= e*e), Schinder & Bludman 



The stationary-state Eddington factor as computed from 
the Boltzmann equation in our program is displayed as a 
function of radius in Fig. The deviation from the rela- 
tivistically correct results of ^chinder fc Bludman (1989) 
is significant. On the other hand, our result is very close to 
the corresponding "weak-field" case, which is not shown in 
Fig. ^ This indicates that general relativistic ray bending, 
which is not correctly taken into account in our calcula- 
tion, has a determining influence on the Eddington factor, 
while the contraction of rods and time dilation — both in- 
cluded in our calculation of the Eddington factor — are 
of minor importance. Indeed, the relativistically correct 
characteristic curves deviate considerably from straight 
rays, as can be seen in [Schinder fc Bludman] (|l989| . Fig. 9). 



(1989) simplified the general relativistic moment equa- 
tions (Eqs. HJEI) to 



(69) 



- — J + — — fR^iJe*) - e*r7(°) 
cBt + R^dR^ uc )-e L, , 

ID ra,,*, s> K ~ J 



9e* 



(70) 



The new coordinates have the advantage that the moment 
equations decouple in energy space (if energy-changing 
neutrino-matter interactions are ignored). Hence, all ra- 
diation quantities depend only parametrically on Eoo- 



The analytical stationary-state solution for the ra- 
diation flux density, which is easily verified^ to read 
i?2e*i/(i?, eoo) = A{eoo), (with A(eoo) e R), is repro- 
duced with an accuracy of 10"''. This, of course, was to be 
expected for the atmosphere with isoenergetic scattering, 
since in the stationary state, the solution H{r) is solely 
determined by the zeroth moment equation (Eq. |69| ) with 
no reference to the (incorrect) Eddington factor. 

More remarkably, the quality of our solution for the 
mean intensity J{r), which is governed in the stationary 
state by Eq. (|70| ) and thus is directly influenced by the 
Eddington factor, appears to be rather good, too. As can 
be seen from Fig. the slope of the general relativistic 
result is reproduced correctly throughout the atmosphere. 
The quantitative agreement is quite satisfactory as well. 
In judging the accuracy of our results one has to keep in 
mind that the model computed here is a rather extreme 
case in two respects: Firstly, the deviations of the metric 
functions e* and F from unity (see Fig. |^) are larger than 
those typically encountered in supernova simulations when 
the neutron star is not going to collapse to a black hole. 
Secondly, as stated in the beginning, the effect of a curved 



spacetime on the radiation field was observed by 3chindei 



in their most general torm (i^qs. p4|, pdj), our implemen- 
tation of the tangent ray scheme for computing the vari- 
able Eddington factors employs a number of approxima- 
tions. In particular, by using straight lines for the tangent 
rays, general relativistic ray bending is not included. The 
current tests therefore serve the purpose of clarifying the 
influence of these approximations on the quality of the 



Different trom the moment equations, which we treat fc Bludman) (|1989D to be much larger for the scattering 



solutions of the moment equations. Schinder fc Bludman 



(1989) found signiflcant differences of the Eddington fac- 
tors for the general relativistic and Newtonian simulations 
of the scattering atmosphere (while in case of the homoge- 
neous sphere the differences were minor). The scattering 
atmospheres therefore seems to be an ideal test case for 
our purposes. 



Following Schinder 



parameters R„ 



= 10 and 



Bludman (1989), we choose the 
The variation of 



= 2. 



the metric functions c and F with radius is depicted in 
Fig. I^a. All other model parameters like boundary condi- 
tions and initial conditions are the same as in the "weak- 
field" -case which is defined by e* = F = 1 (see above). 



atmosphere than in a case dominated by absorption (and 
emission). We expect real situations to be somewhere in 
between these two extremes. 



4.2. Radiative transfer in stationary background media 

In the following section we consider radiative trans- 
fer problems in moving, yet stationary background me- 
dia. This means that in addition to time-independent 
radial profiles of the opacity and emissivity, a time- 
independent velocity field as a function of radius is pre- 
scribed. Stationary-state solutions for the radiation field 
are expected to exist is such cases and have been com- 
puted to high accuracy for some test problems including 
differentially moving atmospheres with relativistic veloci- 
ties. By comparison with fully (special) relativistic calcu- 



^ Note that the more familiar result R^e^"^ H{R) = const, 
only holds for the total (i.e. energy integrated) "flux density" 
H{R) ~ J H{R, e) Ae = e"* / H{R, e^) de^. 
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Metric functions, model GR 
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Fig. 5. Metric functions F (dashed line) and e* 
(solid line) of the background model (Panel a) 
and stationary-state solutions for selected radiation 
quantities ((Panels b, c) versus radius. Panel b shows 
the zeroth moment J (solid line) and first moment 
H (dashed line) of the intensity normalized to the 
non-relativistic solutions. The reference solutions of 
Schinder fc Bludman] ( |198S| ) are given by the sym- 



bols (J: triangles; H: diamonds). Error bars indicate 
the estimated uncertainties of reading off the refer- 
ence solutions from the plots in ^chindcr fc Bludman 



(1989). The Eddington factor vs. radius is shown as 



a solid line in Panel c toget her with the solution of 
S chinder fc Bludman (1989, Fig. 13; diamonds). 
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lations, we are not only able to test our implementation 
of the velocity dependent terms but can also judge the 
quality of the employed C'(?;/c)-approximation of special 
relativistic effects in the equations of radiative transfer. 

Upon introducing a non- vanishing velocity field a par- 
ticular frame has to be specified, where physical quantities 
are measured. In all cases considered in this work, the lat- 
ter is chosen to be the Lagrangian or comoving frame of 
reference. This has to be distinguished from the (numeri- 
cal) coordinate grid that is used to simply label the events 
in spacetime. Although the transformation between differ- 
ent coordinate grids is trivial from an analytical point of 
view (e.g., the simple replacement d/dt -f vd/dr D/Dt 
transforms from Eulerian to Lagrangian coordinates), the 



numerical treatment can be involved (cf. Sect. 3^, where 
our algorithm for computing a formal solution of the ra- 
diative transfer equation in Eulerian coordinates is de- 
scribed). Therefore we present test results obtained with 
two different versions of our radiative-transfer code, one 
which uses Lagrangian coordinates and another one which 
employs an Eulerian coordinate grid. 



Differen 



iaily expanding gray scattering atmospheres: 



Mihalas (198C) presents stationary-state solutions for the 
same type of purely scattering atmospheres considered 
above. In addition to the static case, he investigated differ- 
entially expanding atmospheres with relativistic velocities. 
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Fig. 6. Stationary-state angular mean J (times r^) of the intensity as measured in the comoving frame of reference. 
The abscissa gives the radial position. Results are shown for gray scattering atmospheres that expand according to 
a linear (Model "LVL"; Panel a) and the quadratic (Model "QVL"; Panel b) velocity law. The thin vertical line 
marks the upper boundary of the atmospheres. Different line styles of the curves correspond to different values of 
the parameter /^max which gives the maximum expansion velocity in units of c reached at the surface (r = 11) of the 
model atmospheres. The thin solid curve corresponds to the static case. Reference solutions (symbols) were taken from 



Mihalas (1980, Fig. 2) 



We refer to this paper and study two classes of atmo- 
spheres as test problems, which differ in the functional 
dependence of the expansion velocity {v = c[3) on the 
radius: 



with m e {1, 2} . 



(71) 

The case m = 1 describes a "linear velocity law" , the case 
m = 2 is called "quadratic velocity law" . Different models 
are labeled by the parameter < /3max < 1, which is the 
maximum expansion velocity in units of c reached at the 
surface of the atmosphere (r = rmax)- 

The scattering opacity of the expanding atmospheres 
is the same as the one used for the static cases in Sect. 4.1 
except for the unit of length, which, for the ease of com- 
parison, is adopted from Mihalas ( 19801 ): 



Xit, r, e) = Ks{r) = ar , r^in <r < r,„ax , 



(72) 



= 1/11) in the units of Hummer & 
who set a = 1 (cf . Eq. pq) . For the profile 



with the parameters r^i^ = 1, rmax — H a-nd a — 10.989. 
This yields the same optical depth as the atmosphere with 
'"max = 1 (and r, 

Rybick pl (|l97l| ) _ 
of Eq. (|72|), an optical depth of unity is reached at a ra- 
dius of r w 5.5. Since the opacity does not depend on the 
frequency of the radiation, the atmospheres are referred 
to as "gray" . In this case, the transport solution does not 



depend on e and the 9/9e-terms in the comoving- frame 
Boltzmann equation can be dropped. 

For our simulations we used a numerical grid with 200 
radial points, which were initially distributed logarithmi- 
cally between rmin and rmax- Following Mihalas (1980) as 
closely as possible, we impose the boundary conditions 



-^{i: ^min ; M) 



11 
0, 



< M < 1, 
-1 < < 0, 



(73) 



with R(t) = r^ax + Anax ct. The numerical values in the 
diffusion ansatz at the inner boundary are chosen ac- 



cording to the asymptotic solution of Hummer & Rybicki 
(1971, Eq. 2.22). We started the simulations using the 
stationary-state solution of the corresponding static at- 
mosphere (see Sect. 4.1) as an initial condition (at t = 0). 



The results presented here were obtained by using a 
Lagrangian coordinate grid. Our calculations were termi- 
nated when the radiation field (as measured in the co- 
moving frame of reference) showed no more substantial 
temporal variations at any given radial position. We refer 
to our results as "stationary-state solutions" in this sense. 
We also tested the Eulerian grid version of the code with 
a less extensive set of models. The results were practically 
identical. 

We compare our stationary-state solutions with results 



obtained by Mihalas (198C). In the latter investigation a 



numerical method that is accurate to all orders in v/c 
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Mean intensity, model 1 (m=l, ^q=^0) 



Mean intensity, model 2 (m=1, ^o^^) 
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Fig. 7. Frequency integrated zeroth order angular moment of the comoving frame intensity d^ J{0) ^ function 
of radius for our 0{v/c) stationary-state solutions. The thin vertical line marks the outer boundary of the atmospheres. 
The small arrows point to the positions where an optical depth (which is calculated for ^ < ^o) of unity is reached. 
Reference solutions (symbols) obtained with a fully relativistic code were taken from [Mihalaiij ( |1980| , Figs. 3 and 5), 
where the read-off error is approximately given by the size of the symbols. Different maximum velocities /3max of the 
atmospheres are indicated by different line styles and symbols. In Panel a the results are depicted for Model (1) with 
the opacity step located at = 10 (A = 0.25). Panel b shows results for Model (2) with = 3 (A = 0.2). 



was employed to solve the stationary radiative transfer 
equation. We therefore do not only test the correct nu- 
merical implementation of the 0(v/c) equations (Eqs. ^, 
1^, H), but we can also get a feeling for the quality of the 
C'(u/c)-approximation to the special relativistic transfer 
equation. All characteristic features of the solution, dis- 
cussed in detail by Mihalas| ( 1980| ), are reproduced cor- 
rectly by our implementation of the 0{v/c) equations. It 
is remarkable how accurately the fully relativistic solu- 
tions are reproduced (Fig. ^). The differences are < 10%, 
even for /3max = 0.5. 



DifFeren :ially expanding, nongray, isothermal atmospheres: 



Mihalas (1980) computed stationary-state solutions also 
for relativistically expanding nongray atmospheres, i.e., 
the opacity and emissivity depended explicitly on the ra- 
diation frequency. Therefore the full frequency dependence 
of the comoving frame Boltzmann equation had to be re- 
tained. The atmospheres were assumed to be stationary, to 
emit a thermal continuum of radiation, and to be isother- 
mal, i.e., no radial or temporal variations of the source 
function were present. The radiation-matter interaction 
was solely via absorption and emission. 

Introducing the dimensionlcss "frequency" ^ — e/T 
(the temperature T is measured in units of the Boltzmann 



constant fee ) , the emissivity reads 



e 



exp(0 - 1 



(74) 



The absorption opacity Ka was given a step-like functional 
dependence on frequency, characterized by the frequency 
^0 where the step is located and some characteristic width 
A of the function s(^) := exp[— (^ — ^o)^/A^], which me- 
diates a continuous variation of Ka across the step: 



«a(7-, ?) 



XiM • 5(0 +X2(r) -(1-5(0), fore<eo, 
Xi{r), for^>^o, 

(75) 

with xi(r) := lOar ^, x^i"!") '■= ctr ^, and the value 
a := 10.9989. Thus, the optical depth at a given radius 
is roughly 10 times smaller for "low-frequency" radiation 
{£, < S,o) than for radiation with "high" frequency (^ > ^o)- 

Following Mihalas (1980), we consider the two exam- 
ples (1) Co = 10, A = 0.25, and (2) = 3, A = 0.2. In 
the former case the frequency ^o, where the opacity jump 
is located, is considerably larger than the frequency of the 
maximum of the Planck function (^max ~ 2.82), whereas 
it is close to the maximum in the latter case. 

The radial extent of the atmospheres was chosen to 
be < r < 11 in both cases. For all models the linear 
law (m = 1, cf. Eq. ^l]) for the expansion velocity as a 
function of radius was used. The results presented below 
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Fig. 8. Spectral distributions of the angular moment J of the comoving frame intensity for our stationary-state 
solutions of the 0{v/c) radiative transfer equation. Spectra are displayed for a radial position corresponding to an 
optical depth (which is calculated for ^ < ^o) of unity and for the surface of the atmosphere. For comparison, 
the equilibrium intensity (Planck function) at an optical depth of unity is plotted as a thin line. Different maximum 
velocities (/3max) of the expanding atmospheres are indicated by different line styles. Panel a shows results for Model (1), 
where the opacity step (vertical arrows in the plots) is located at = 10 (A = 0.25). Panel b shows Model (2) with 



^0 = 3 (A = 0.2). Our results may be compared with Mihalae (1980, Figs. 4 and 6) 



were obtained using an Eulerian coordinate grid. Our code 
version with a Lagrangian grid yields very similar results. 



We have used 200 logarithmically distributed radial 
zones for the radial range 1 < r < 11, supplemented by 
a few additional zones between < r < 1. Boundary 
conditions were chosen as usual: Spherical symmetry at 
the coordinate center requires I(t, r = 0, /i = 1) = X(t, r — 
0,fj, — —1), and demanding that no radiation is entering 
through the surface of the atmosphere translates to the 
condition I{t, r = 11, n < 0) — 0. 

The energy grid consisted of 21 bins (very similar re- 
sults were obtained using 8 and 12 bins instead) of equal 
width covering the range < ^ < 12 for Model (1) and 
< ^ < 8 for Model (2), respectively. Details about the 
stationary-state solutions and their physical interpreta- 

pV). In our test cal- 



tion can be found in Mihalas ( 1980 



culations we were able to reproduce all qualitative trends 
discussed there. The quantitative agreement of the zeroth 
order angular moment J of the specific intensity is within 
a few percent (see Fig. ^), even for the case of /3max — 0.5. 
This, as previously stated, does not only demonstrate the 
accuracy of our numerical implementation but even more 
importantly also justifies the physical approximations cm- 
ployed in the underlying finite difference equations. Note 
that in applications of our code to supernova calculations, 
one expects velocities that are not in excess of ~ 0.3c. 



A rigorous comparison of the spectral distribution of 
the angular moment J as calculated by our radiative trans- 
fer code (cf. Fig. ||) with the results of Mihalas (198C) is 
difficult, since he shows only results for the static case 
and the rather extreme case /3max = 0.8. The latter atmo- 
sphere obviously cannot be modeled reasonably well with 
our 0{f3)-code. Nevertheless, by inspection of Figs. 3 and 
4 of Mihalas (198C) one can infer that the shapes of our 
spectra as displayed in Fig. ^ exhibit the correct qualita- 
tive dependence on the expansion velocity. 



Frame dependence? As an important consistency check 
we verified that, to order 0{v/c), the numerical solution 
of our implementation of the comoving frame equations 
shows the correct transformation properties of the stress- 
energy tensor of radiation (see Eq.p^) and no unphysical 
artifacts are caused by the choice of a moving reference 
frame (cf. the critique of Lowrie et al. 1999). 



For this purpose results are produced for neutrino 
transport in a thermally and hydrodynamically frozen 



post-bounce model of a supernova calculation by Brucnn 
(1993). This model can be characterized as follows: 



Neutrinos are emitted from a hot and dense hydrostatic 
inner core with a radius of r w 110 km. This inner core is 
surrounded by a supersonically infalling outer core, which 
is effectively transparent to the radiation. Velocities in this 
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Fig. 9. Comparison of the stationary-state solutions of the comoving frame transfer equation (bold lines) with the 
solutions after transformation to the Eulerian frame (thin lines) and a reference calculation for a static stellar back- 
ground {v = 0; dashed lines). Panel a shows the total energy density of neutrinos (times c), Panel b displays the 
total energy flux density (times r^). 



outer atmosphere range from v{r 
v{r = 250 km) « -0.08c. 



140 km) « -0.12c to 



The stress-energy tensor as derived from the 
stationary-state solution of the comoving frame transport 
equation (Eq. |6|) is transformed to the Eulerian frame 
(i.e. the inertial frame in which the center of the star is at 
rest) and can then be compared with the result of an in- 
dependent calculation which solves the transport equation 
directly in the Eulerian frame. 

Figure ^ shows the components E — Att/c de J and 
F = An Jp°° de H of the stress-energy tensor (scaled by r^c 
and r^, respectively). Due to the large inwardly directed 
velocities in the outer atmosphere, the neutrinos emitted 
by the inner core are blueshifted for observers which arc 
locally comoving with the fluid flow (see the bold lines 
in Fig. ^). When the stress-energy tensor is transformed 
to the Eulerian frame (E q.[lO| ) , this effect obviously disap- 
pears (thin lines in Fig. |9|). Comparison with the stress- 
energy tensor obtained by solving the transport equation 
directly in the Eulerian frame (dashed lines in Fig. H) re- 
veals good agreement for both the energy density E and 
the energy flux density F in the rapidly moving atmo- 
sphere. This demonstrates that no unphysical frame de- 
pendence is present in the calculations. The fluctuations 
of the transformed energy flux density near the center are 
caused by the term /3{J + K) in the expression for i?^"' 
in Eq. (p^). Because J (and K) are large compared with 
H in the dense interior, even small velocities lead to a 
significant "advective contribution" to H^^"^. 



Note that an Eulerian-frame calculation in general re- 
quires complicated velocity-dependent transformations of 
the source terms on the rhs. of the transport equation 
(see Mihalas & Mihalas 1984). This, however, is unnec- 
essary in the specific situation of the discussed model: 
Because the region where most of the neutrinos are pro- 
duced moves with low velocities, the source terms there 
can be evaluated in the rest frame of the fiuid. In the 
outer atmosphere on the other hand, where the large ve- 
locities would require a careful transformation, the source 
terms nearly vanish. This allowed us to simply drop the 
velocity-dependent terms on the Ihs. of Eq. (^) in order to 
perform a calculation with results that are in agreement 
with the Eulerian frame solution in the low-density part 
of the star. 

4.3. Core collapse and supernova simulations 
4.3.1. Newtonian gravity 

Our new neutrino hydrodynamics code has already been 
applied to dynamical supernova simulations (Rampp & 



Janka 2000). Here we report some results of a Newtonian 



calculation of the collapse phase of a stellar iron core with 
mass M-pc — 1-28 Mq (plus the innermost 0.1 Mq of the 
silicon shell). It is the core of a 15 M© progenitor star 
(Model "sl5s7b2", |Woosley||1999|; |Heger| |2000|) . We com- 



pare the results to a calculation published by Bruenn & 



Mezzacappa (1997). The input physics, in particular the 



stellar model and the high-density equation of state, were 



adopted from model "B" of Bruenn & Mezzacappa (1997), 
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Fig. 10. Profiles of the total lepton fraction Yiep (dash- 
dotted line) and the electron neutrino fraction Y,y^ (solid 
line) as a function of the enclosed mass at the time when 

For comparison 
Fig. 5a) are 



the central d ensity reaches 10"'^'^ Kcm ^. 



the results of Bruenn fc Mezzacappa ( 1997 



plotted with diamonds. 

with the exception that we do not include neutrinos other 
than Vc- This, however, does not make a difference dur- 
ing the collapse phase where the electron degeneracy is so 
high that only Vc can be produced in significant numbers. 
The hydrodynamics was solved on a grid with 400 radial 
zones out to 20 000 km, which were moved with the in- 
falling matter during core collapse. For the transport we 
used an Eulerian grid with 213 geometrically spaced radial 
zones, 233 tangent rays and 21 energy bins geometrically 
distributed between and 380 MeV the center of the first 
zone being at 2 MeV. 

The total lepton fraction Yicp :^ + Y^^ and the en- 
tropy per baryon at the center of the core as a function of 
the central den sity are displayed in Fig . [Tl|. T he agreement 
with results of Bruenn fc Mezzacappa ( 1997 ) is excellent. 
Also the total lepton fraction as a function of enclosed 
mass matches perfectly (Fi g. [l0|). Defining the sho c k for - 
mation point according to Bruenn & Mezzacappa (1997) 
as the mass coordinate where the entropy per baryon af- 
ter core bounce first reaches a value of s = 3 fee we find 
Ms w O.62M0. This is about O.O3Af0 or 5% larger than 



the value given by Bruenn & Mezzacappa (1997, Table V). 

Judging our results, one should, of course, take into 
account that Bruenn & Mezzacappa (1997) employed 



a multi-group flux-limited diffusion (MGFLD) approx- 
imation for treating the neutrino transport, whereas 



our code solves the Boltzmann equation. Mezzacappa & 
Brueni| ( |1993b| J^ performed a detailed comparison of 
core-collapse simulations with MGFLD and with their 



Boltzmann solver based on the SN-method. For the quan- 
tities presented here they also found good agreement of 
their results throughout the inner core. 

In the outer regions of the collapsing core, the situation 
is somewhat different: Mezzacappa Sz Bruenn] ( |1993b| ,|c|) 



demonstrated that a number of quantities reveal signif- 
icant deviations between MGFLD and the Boltzmann 
transport (SN-method). In particular, they obtained an 
electron neutrino fraction which was by up to 20% larger 
in the Boltzmann run. For the total lepton fraction, how- 
ever, the difference was only 1%. Our results (Fig. |l^ con- 
firm the agreement of Boltzmann and MGFLD results for 
the latter quantity. A more detailed comparison between 
the variable Eddington factor method and the SN-method, 
however, is desirable and currently underway. 

Spectral resolution: A few quantities exhibit spurious os- 
cillations as a function of the radial (mass-) coordinate and 
time in our core-collapse simulations. Figure 12 a, shows 
such features for the electron neutrino fraction in the 
mass shells < M < 0.6 Mq when the central density 
is pc — 10^* gcm"'^. Similar structures are visible in the 
profile of the entropy per baryon and in the neutrino lu- 
minosity as a function of radius (or mass) once the central 
density becomes larger than pc ~ 10"'^^ gcm~'^. They can 
also be seen in the central entropy as a function of the cen- 
t ral de nsity or time (cf. Fig. pTpp). Mezzacappa & Bruenn 
( 1993c ) found a similar effect in their core-collapse calcu- 
lations. They identified the finite spectral resolution of the 
degenerate Fermi distribution of t he electron neutrinos to 
be the origin of these oscillations ( Mezzacappa fc Bruenn 
1993c| , §3.3). 



For comparison, we have computed a collapse model 
with a significantly improved spectral resolution: Instead 
of 21 bins we used 81 energy zones to describe the energy 
spectrum between and 380 MeV. A result of this simu- 
lation is displayed in Fig. The neutrino fraction (and 
related quantities) are now perfectly smooth functions of 
the mass coordinate. 

A reasonable compromise between accuracy and com- 
putational load could be achieved by using approximately 
20-30 (geometrically spaced) energy bins. This seems to 
be sufficient to adequately treat the sharp Fermi surface 
of the degenerate distribution of electron neutrinos that 
builds up during collapse. Even for less energy bins fluc- 
tuating quantities follow the correct evolutionary trends 
(Fig. [l^a). Therefore we conclude with the remark that 
the presence of oscillations in some transport quantities 
may be considered as undesirable. The overall evolution 
of the collapsing core, its deleptonization or the shock for- 
mation radius, however, were found not to be sensitive to 
such "noise". 

Also during the post-bounce evolution (for which the 
transport of i>c was included), we determined no signiflcant 
differences between the dynamical evolution of a model 
with 17 energy bins and a model computed with 27 en- 
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Fig. 11. Central lepton fraction Y|op :— Yc + Y^,^ (Panel a) and central entropy per baryon (b) as functions of the 
central density for our Newtonian core-collapse simulation of a 15 M© star (solid lines), compared with the results 
(diamonds) obtained by Brucnn & Mezzacappa (1997, Figs. 4a, 8a). 



ergy bins, using an otherwise identical setup: At a time 
of 100 ms after bounce, the structure of the models, mea- 
sured by the radial positions of fluid elements with the 
same mass coordinate, agreed within a relative error of 
< 1%. This is corroborated by the fact that for given snap- 
shots from the dynamical evolution we found the neutrino 
source terms Qe and also Qn to be practically identical 
in both calculations. 

Radial resolution and size of the numerical time step: 

Varying the number of radial zones for the hydrodynamic 
and transport parts of the code between 100 and 200 and 
using a comoving hydrodynamic grid during the phase of 
stellar core collapse, we found no significant changes of the 
results. The same holds true for the post-bounce evolution 
where we switch to a fixed, Eulerian grid. The only excep- 
tion to this insensitivity is the central density at bounce, 
which shows variations of up to 10%, depending on the 
size of the innermost zones of the hydrodynamic grid. 

The choice of the radial grid, however, also influences 
the size of the numerical time step. According to our expe- 
rience the latter has to be chosen very carefully. We give 
two examples here: 

(a) While neutrinos escape unhindered initially, the col- 
lapsing stellar core becomes opaque when the density 
exceeds ~ 10^^ gcm^'^, and "neutrino trapping" sets 
in. Since neutrinos are then carried along with the stel- 
lar plasma, the subsequent, accelerating phase of the 
infall is adiabatic (the sum of medium and neutrino 
entropy per baryon is constant: s + = const.) and 
the total lepton fraction of a fluid element is conserved. 



Fulfilling the corresponding conditions, DYJep/Dt = 
and D(s -j- s^)/T)t = 0, numerically on an Eulerian 
grid requires a sufficiently accurate treatment of the 
advection of energy and leptons. We found that treat- 
ing the advection terms in the neutrino transport at 
first order accuracy in time (setting ^ = 1 in Eq. p5| ) is 
not sufficient and causes deviations from i^ep = const, 
and s + Si, = const, at a level of > 10%. This prob- 
lem can be diminished when the transport time step 
is reduced significantly below the limits mentioned in 



Sect. 3.6. The time step reduction and the correspond- 



ing increase in computer time, however, can be avoided 
by handling the advection "second-order" accurate in 
time, chosing C « 0.5 in Eq. (|2^). As visible from 
Fig. In], this yields very good results (note that the en- 
tropy of neutrinos, Si/, is small compared to the plasma 
entropy s, and essentially constant after neutrino trap- 
ping). 

(b) Another problem occurred during the long-time evolu- 
tion of the forming neutron star after core bounce. For 
too large a value of the hydrodynamic time step, we 
discovered an artificial inflation of the mass shells that 
are located in the region of a steep density gradient 
below the neutrinospherc. This effect was found not 
to be caused by insufficient energy resolution of the 
neutrino transport nor does it directly depend on the 
employed number of radial zones. Varying the latter 
by a factor of two and testing 17 and 25 energy bins 
in the interval between and 380 MeV, we found that 
the agreement of our transport results was better than 
1% and thus not worse than the uncertainties due to 



30 



M. Rampp and H.-Th. Janka: Radiation hydrodynamics with neutrinos 
Neutrino fraction, model WS18_lr Neutrino fraction, model WS18_hr 



0.08 





Fig. 12. Electron neutrino fraction Y^^ (bold solid line, scale given by the axis on the left side of the plots) and 
electron neutrino equilibrium chemical potential /i^'^ ("Fermi surface", dash-dotted line, scale given by the axis on 
the right side) as functions of enclosed mass at the time when the central density has reached lO^"* gcm~^(The plots 



show results of collapse simulations for a ISM© blue supergiant progenitor star; Woosley et al. 1997). The spectral 



resolution of the energy grid is indicated by the spacing of the short horizontal lines (right ordinate), which correspond 
to the boundaries of the energy grid. Panel a displays results obtained with "low" spectral resolution (21 energy bins 
distributed geometrically between and 380 MeV), Panel b shows the run with "high" spectral resolution (81 energy 
bins covering the same spectral range) . For better comparison YJy^ of the high- resolution run is drawn as a thin solid 
line also in Panel a. 



the interpolation of the stellar background. Since the 
post-bounce evolution is usually calculated with the 
same grids for the hydrodynamics and the transport, 
possible errors associated with the mapping of quan- 
tities between the grids can also be excluded as the 
source of the problem. Moreover, even if different ra- 



dial grids are used as in the model published by Rampp 



& Janka (200C), an expansion of the outer layers of the 
forming neutron star is not necessarily a consequence. 
The quality of our calculations, however, turned out 
to be dependent on the size of the time step of the hy- 
drodynamics code. Our tests showed that a reduction 
of the time step was necessary to prevent the artifi- 
cial expansion of the quasi-hydrostatic structure. We 
were so far unable to exactly locate the origin of this 
numerical problem, but presume that it might be con- 
nected with the accuracy of the implementation of the 
neutrino momentum term in the hydrodynamics code. 
Good results, however, were obtained, when the length 
of the time step was guided by the typical time scale 
set by the neutrino momentum transfer to the stellar 
medium, Atnyd — pv/Qu (cf- Eq. ||). 



4.3.2. Approximate relativistic core collapse 

We also tried to assess the quality of our approximate 
treatment of general relativity (GR) in the equations of 
hydrodynamics and neutrino transport (cf. Sect. 3.7). 



For this purpose we performed core-collapse simulations 
for the 15 Mq progenitor star "sl5s7b2" and compared 
characteristic quantities with the information given for a 
MGFLD simulation of the same star in a recent paper by 



Licbcndorfcr ct al 



Brucnn ct al. (2001 ), and for a Boltzmann simulation by 



( |2001aD . 

For the hydrodynamics we used 400 radial zones out 
to 10 000 km, which were moved with the infalling mat- 
ter during core collapse. The transport was solved on an 
Eulerian grid with 200 geometrically spaced radial zones, 
220 tangent rays and 17 energy bins geometrically dis- 
tributed between and 380 MeV the center of the first 
zone being at 2 MeV. 

At bounce, which occurs at t^^ = 177.7 ms (i^"'"* = 
196.1 ms) after the start of the calculation, the central 
density reaches a maximum of = 4.8 x 10^"* gcm^'^ 
(pNowt ^ 3 35 ^ iQii gcni-3) in our GR model (for com- 
parison, the corresponding values for the Newtonian runs 
are given in brackets). The hydrodynamic shock forms at a 

0.51 Mq {M^''^^ = 0.62 Mq). 



mass coordinate of Mg^ 
For the same stellar model, iBruenn et alJ (pOOll) found a 
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bounce density of = 4.23 x 10^^ gcm^'^ for the gen- 
eral relativistic simulation (p^''** = 3.20 x 10^'' gcm^'^). 
The corresponding time of the bounce was = 182.9 ms 
(^Newt ^ 201.3 ms). piebendorfer et al.| ( p001a|) located the 
shock formation at a mass coordinate of M^^ = 0.53 Mq 
(M^™* = 0.65 Mq). 

Since the collapse time is determined by the slow ini- 
tial phase of the contraction, which is sensitive to the 
initial conditions and the details of the numerical treat- 
ment at the start of the calculation (Liebendorfer, per- 
sonal communication), it may be better to use the dif- 
ference between the GR and the Newtonian simulations, 
instead of comparing absolute times at bounce. Our value 
18.4 ms is very close to the result of 

j-Nowt 



of 



Bruenn et al.| ( |200lD : t^"""* - = 21.6 ms. A part of the 
remaining difference may be attributed to the fact that 
Bruenn et al. (2001) employed a MGFLD approximation 
of the neutrino transport with some smaller differences 
also in the input physics (e.g., Bruenn et al, 2001 as well 



as Liebendorfer et al, 2001a did not take into account ion- 



ion correlations for the coherent scattering, but we did). 

At the level of comparison which is possible on grounds 
of published numbers we conclude that the agreement be- 
tween our approximate treatment of GR and the relativis- 
tic core collapse seems to be reasonably good. 

5. Summary and conclusions 

We have presented a detailed description of the numerical 
implementation of our new neutrino transport code and its 
coupling to a hydrodynamics program which allows sim- 
ulations in spherical symmetry as well as in two or three 
dimensions. The transport code solves the energy and time 
dependent Boltzmann equation by a variable Eddington 
factor technique. To this end a model Boltzmann equa- 
tion is discretized along tangent rays and integrated for 
determining the variable Eddington factors which provide 
the closure relations for the moment equations of neutrino 
energy and momentum. The latter yield updated values of 
the neutrino energy density and neutrino flux, which are 
fed back into the Boltzmann equation to handle the colli- 
sion integral on its right hand side. 

The system of Boltzmann equation and moment equa- 
tions together with the operator-splitted terms of lepton 
number and energy exchange with the stellar background 
(which influence the evolution of the thermodynamical 
quantities and the composition of the stellar medium and 
thus the neutrino-matter interactions) are iterated to con- 
vergence. Lepton number conservation is ensured by solv- 
ing additional moment equations for neutrino number den- 
sity and number flux. 

The integration is performed with implicit time step- 
ping, which avoids rigorous limitations by the CFL condi- 
tion and ensures that equilibrium between neutrinos and 
stellar medium can be established accurately and without 
oscillations despite of stiff neutrino absorption and emis- 
sion terms. The coupling of energy bins due to Doppler 
shifts and energy changing scattering processes is imple- 



mented in a fully implicit way. When coupling the trans- 
port part to an explicit hydrodynamics program, as done 
in this work, computational efficiency requires to have the 
option of using different time step lengths for transport 
and hydrodynamics. In addition, we found it advantageous 
to have implemented the possibility of choosing different 
radial grids for both components of the program, and to 
switch between Lagrangian and Eulerian coordinates de- 
pendent on the considered physical problem (although we 
always measure physical quantities of the transport in the 
comoving frame of reference) . 

We have suggested an approximation for applying the 
code to multi-dimensional supernova simulations which 
can account for the fact that hydrodynamically unstable 
layers develop in the collapsed stellar core. The variable 
Eddington factor method offers the advantage that this 
can be done numerically rather efficiently (concerning pro- 
gramming effort as well as computational load) . Although 
the neutrino moment equations for the radial transport 
are solved independently in all angular zones of the spher- 
ical grid of the multi-dimensional hydrodynamical simu- 
lation, the set of moment equations is closed by a variable 
Eddington factor which is computed only once for an an- 
gularly averaged stellar background. 

This approximation saves a significant amount of com- 
puter time compared to truely multi-dimensional trans- 
port, and yields a code structure which can easily be im- 
plemented on parallel computers. The treatment is con- 
strained to radial transport and thus neglects lateral prop- 
agation of neutrinos. Its applicability is therefore limited 
to physical problems where no global asphericities occur. 
It must also be considered just as a first, approximate step 
towards a really multi-dimensional transport in convective 
layers inside the neutrinosphere. The approach, however, 
should be perfectly suitable to treat anisotropics associ- 
ated with convective overturn in the neutrino-heated re- 
gion between the neutron star and the supernova shock. In 
the latter region neutrinos are only loosely coupled to the 
stellar medium (the optical depth is typically only 0.05- 
0.2). Therefore the energy and lepton number exchange 
between neutrinos and the stellar medium depends on the 
presence of anisotropics and inhomogeneities, whereas the 
transport of neutrinos is essentially unaffected by non- 
radial motions of the stellar plasma. 

In preliminary multi-dimensional simulations of trans- 
port in convecting neutron stars we have convinced our- 
selves that this method guarantees good numerical con- 
vergence, because the variable Eddington factor as a nor- 
malized quantity depends only weakly on lateral varia- 
tions in the stellar medium. The method ensures global 
energy conservation and enables local thermodynamical 
equilibrium between stellar medium and neutrinos to be 
established when the optical depth is sufficiently high. We 
therefore think that the proposed approximation is prac- 
ticable for neutrino transport in multi-dimensional super- 
nova simulations before fully multi-dimensional transport 
becomes feasible. The latter is certainly a major challenge 
for the years to come. 



32 



M. Rampp and H.-Th. Janka: Radiation hydrodynamics with neutrinos 



Our transport code exists in a 0{v/c) version for non- 
relativistically moving media, and in a general relativis- 
tic version. We have not yet coupled it to a hydrody- 
namics program in full general relativity. Instead, we per- 
formed calculations of approximate relativistic core col- 
lapse, where corrections to the Newtonian gravitational 
potential were included and only the redshift and time di- 
lation effects were retained in the transport equations (this 
means that space-time is considered to be flat). The results 
for stellar core collapse compared with published fully rel- 
ativistic calculations (with Boltzmann neutrino transport 
by Liebendorfer et al. 2001 and with multi-group flux- 
limited diffusion by Bruenn et al. 2001) are encouraging 
and suggest that this approximation works reasonably well 
and accounts for the most important effects of relativistic 
gravity as long as one does not get very close to the limit 
of black hole formation. 

We have presented results for a variety of idealized, 
partly analytically solvable test calculations (in spher- 
ical symmetry) which demonstrate the numerical effi- 
ciency and accuracy of our neutrino transport code. The 
"neutrino radiation hydrodynamics" was verified by core- 
collapse and post-bounce simulations for cases where pub- 
lished results were available and could reasonably well 
serve for a comparison (a direct and more detailed com- 
parion with Boltzmann calculations using the Sn method 
by M. Liebendorfer and A. Mezzacappa is in progress). 

Tests for a number of static model atmospheres showed 
that the treatment of the angular dependence of the neu- 
trino transport and phase space distribution can be han- 
dled accurately by the variable Eddington factor method. 
This holds for moderately strong general relativistic prob- 
lems, too, even if ray bending effects are neglected in the 
determination of the variable Eddington factor. The nu- 
merical quality of the handling of the energy dependence 
of the transport was also checked by considering back- 
ground models with stationary fluid flow. We included a 
model with mildly relativistic motion (up to v/c = 0.5) 
and found that the code produces accurate and physically 
consistent results although we employ an 0{v / c) approx- 
imation to the special relativistic, comoving frame radia- 
tive transfer equation. Also the omission of some velocity- 
dependent terms in the model Boltzmann equation for de- 
termining the variable Eddington factor turned out not to 
be harmful in this respect. Closing the radiation moment 
equations by a variable Eddington factor seems to be re- 
markably robust against approximations in the treatment 
of the Boltzmann equation. 

The tests have also demonstrated that our code per- 
forms very efficiently. Only a few iterations (typically be- 
tween 3 and 7) are needed for obtaining a converged so- 
lution of the system of Boltzmann equation and moment 
equations even in cases where scattering rates dominate 
neutrino absorption. The calculation of the formal solu- 
tion of the Boltzmann equation, from which the variable 
Eddington factor is derived, turned out to consume only 
10-20% of the computer time. The major computational 
load comes from the inversion of the set of moment equa- 



tions. The latter have a reduced dimensionality relative 
to the Boltzmann equation, because the dependence on 
the angular direction of the radiation propagation was re- 
moved by carrying out the angular integration. This ad- 
vantage, however, has to be payed for by the iteration 
procedure with the Boltzmann equation. 

Using Boltzmann solvers for the neutrino transport 
means a new level of sophistication in hydrodynamical 
simulations of supernova core collapse and post-bounce 
evolution. It permits a technical handling of the trans- 
port which for the first time is more accurate than the 
standard treatment of neutrino-matter interactions. The 
latter includes a number of approximations and simplifi- 
cations which should be replaced by a better description, 
for example the detailed reaction kinematics of neutrino- 
nuclcon interactions, phase space blocking of nucleons, ef- 
fects due to weak magnetism in neutrino-nucleon reac- 
tions, or nucleon correlations in the dense medium of the 



neutron star (Reddy et al. 1998 , 1999 ; Burrows 
1999| ; |Horowitz| ^0"02|) ! 



Sawyer 



Boltzmann calculations will not only remove impon- 
derabilities associated with the use of approximate meth- 
ods for the neutrino transport in hydrodynamical super- 
nova models. They will also allow for much more reliable 
predictions of the temporal, spectral and flavour prop- 
erties of the neutrino signal from supernovae and newly 
formed neutron stars. This is an indispensable require- 
ment for the interpretation of a future measurement of 
neutrinos from a Galactic supernova. 
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Appendix A: Neutrino opacities 

In this appendix we shall describe the neutrino-matter 
interactions that are included in the current version of 
our neutrino transport code for supernova simulations. We 
shall focus on aspects which are specific and important for 
their numerical handling and will present final rate expres- 
sions in the form used in our code and with a consistent 
notation. 
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Reaction Rate described in Reference 


i^e* ^ i/e* Sects. 
vA ^ vA Sects. 

!/N Sects. 
VcTi ^ e~p Sects. 
VeP ^ e^n Sects. 
VeA' ^ e~A Sects. 
uu ^ e~e^ Sects. 
wNN ^ NN Sects. 


A. 1.2, 


A. 2.4 




MczzacapDa ,t BruemJ (|19931J): pcrnohorskvl (11994) 


A. 1.2, 


A. 2. 5 




Horowitz 


(L997 


; Bruenn & Mezzacappa 


(19971) 


A. 1.2, 


A. 2.] 


Brucnn 


( 
( 
( 
( 
( 


1985 


); 
); 
); 
); 
); 


Mezzacappa & Brucnn 


( 
( 
( 
( 


1993c 


) 
) 
) 
) 


A.1.1, 


A. 2.1 


Brucnn 


1985 


Mezzacappa <t Brucnn 


1993c 


A.1.1, 


A. 2.1 


Brucnn 


1985 


Mezzacappa & Brucnn 


1993c 


A.1.1, 


A. 2. 2 


Bruenn 


1985 


Mezzacappa & Bruenn 


1993c 


A. 1.3, 


A.2.E 


Brucnn 


1985 


Pons et al (|199^) 








A. 1.3, 


A.2.e 


Hanncstad & Raffcit (1998) 





Table A.l. Overview of all neutrino-matter interactions currently implemented in the code. For each process we list 
the sections where we summarize fundamental aspects of the calculation of the corresponding rate and give details 
of its numerical implementation. The references point to papers where more information can be found about the 
approximations employed in the rate calculations. In the first column the symbol v represents any of the neutrinos 
i^o, i^cj i'/ii i^fj.! t^rj J'tj the symbols e~, e"*", n, p and A denote electrons, positrons, free neutrons and protons, and heavy 
nuclei, respectively, the symbol N means n or p. 



A.l. Basic considerations 



A.1.1. Neutrino absorption and emission 



In the following we discuss the different neutrino-matter 
interaction processes which contribute to the source term 
("collision integral") on the rhs. of the Boltzmann equa- 
tion: 





cat or r 



{1} 



(A.l) 



The sum on the rhs. of this equation runs over all consid- 
ered interaction processes that can change the distribution 
function of a particular type of neutrino. For simplicity we 
have written down the Boltzmann equation for static me- 
dia, here. In the general case of a moving stellar fluid, 
velocity-dependent terms have to appear on the Ihs. of 
Eq. (A.l) when the neutrino quantities and the neutrino- 
matter interaction rates are evaluated in the local rest 
frame of the fluid. The phase-space distribution function 
/ is related to the specific intensity X used in the main 
body of this paper by X = c/{2Trhc)^ ■ f. Accordingly, 
the source term on the rhs. of the Boltzmann equation, 
Eq. (|6|), and the corresponding moment equations for neu- 
trino number and neutrino energy or momentum (Eqs. 
|3l| , 1^, ^) can be calculated from Bj as 



(7(*^)(e) 
CC^He) 



c 
c 



e3.^i?,(e,^), 
{1} 

{1} 

{1} 



where B'^^\e) := 1/2 J^^ dfi ii'' Bi{e, fi). Here and in the 
following, the dependence of quantities on the space- 
time coordinates (t, r) is suppressed in the notation. 
Throughout Appendix ^ the temperature T is measured 
in units of energy. 



The rate of change (modulo a factor of 1/c) of the neu- 
trino distribution function due to absorption and emission 
processes is given by (sec Bruenn 1985| ) 



BAE{e,fi) = j{e)[l - f{e,^^)] - f{e,fi)/X{e) , (A.3) 

where j denotes the emissivity and A is the mean free path 
for neutrino absorption. The factor (1 — /) accounts for 
fermion phase space blocking effects of neutrinos. Using 
the Kirchhoff-Planck relation ( "detailed balance" ), and in- 
troducing the absorption opacity corrected for stimulated 
absorption. 



1 



1 



1 



1 - ' A A ' 



Eq. (A.3) can be rewritten as 



BAB{e,^^)^n:{e)[^{e)~f{e,^l)]. 



(A.4) 



(A.5) 



In Eq. (As) it becomes evident that the source term drives 
the neutrino distribution function /(e, fj.) towards its equi- 
librium value /°'^(e) = (1 -I- exp [(e — /i°'^)/r])~-^, where 
fj,°^ denotes the chemical potential for neutrinos in ther- 
modynamic equilibrium with the stellar medium. In equi- 
librium, the source term vanishes in accordance with the 
requirement of detailed balance. 



(A. 2) A. 1.2. Scattering 

Reduction of the collision integral: The rate of change of the 
neutrino distribution function due to scattering of neu- 
trinos off some target particles is given by the collision 
integral (cf. |Cernohorsky||l994 Eqs. 2.1, 2.3) 



Bsiq) 



(2^ ■ 
[/:(l-/.)ff"(q,q') 



/.(l-/:)i?°"'(q,q')] 



(A.6) 
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where f[, depends on q' and the scattering kernels i?"^ 
and are defined as the following phase space inte- 

grals over products of the transition rate TZ and the Fermi 
distribution functions _Fx of the target particles: 



i?°"*(q,q') = 



C • {2TThc)^ 

(2^(2^ 
2 

c • {2T:hc)^ ' 
d^p' d^p 

(27r)3 (27r)3 



{l-FT)FJ,TZ{p',q';p,q), 



{l-F^)FTTZ{p,q;p',q'). 



(A.7) 



The vectors q and q' are the momenta of the ingoing and 
outgoing neutrinos, respectively, and p and p' those of 
their scattering targets. 

The scattering kernels are usually given as functions of 
the energies e and e' of the ingoing and outgoing neutrino 

and the cosine uj = j.ifi' + — — /i'^) cos(iy9 — (p') 

of the scattering angle, where (/i, tp) and (/x', ip') are the 
corresponding momentum space coordinates. Expanding 
the scattering kernels (i?"Vout i/(27rc)3 • i^in/o"*^ -^^ 
Legendre series 



i?'"/°"'(e,e',w) = 



21 + 1 



1=0 



with the Legendre coefficients 



, in/out/ /X 



dtjPi(u;)i?'"/°"*(e,e',t^), (A.9) 



and applying the addition theorem 
Legendre polynomials Pi{^^) — Pi{^)Pi{fJ.') + 
2ELi fe^PrMJ^"^(A^; )cos[m(^ - (^')] (e.g., 
Bronstein fc Semendjajew 1991 , P™ are the associ- 
ated Legendre polynomials), the integral over tp' in the 



collision integral can be performed analytically ( Yuch fc 



Buchlc| |1977| ; [Schinder fc Shapiro| |l982| ). For use m our 
Boltzmann transport code, it is necessary to truncate 
the Legendre expansion at some level. Note that the 
ortho gonality relation f^^ duj Pi {lu)P ii (uj) — 2lTT'^"' 
(e.g., Bronstein & Semendjajew 1991) implies that any 



truncation of the Legendre series i?'"/°"'(e, e', w) still 
gives the exact integral value J^^ duj R^'^^°^^{e,e' ,uj) 
independent of the level of truncation /max > 0. 



The collision integral and its first two angular moments 
finally readn: 



de' e' 



{< 



1 - 



Bs(e,M) = 27r 

CXD 

/.) . ^(2/ + l)PKM)C(e,e'W(e') 



^-0 



1=0 



(A.IO) 



1 / /2 

de 6 



CX3 

{j2{2l + l){Sio-Li{e)W{e,e')L[{e') 

oo 

J2{2l + l)Li{e^r\^,^'Mo-me'))} 



(A.ll) 



Bl}\e)^2Tr / de' e' 



{ J2[Sn - (/ + l)P;+i(e) " lUMe)W{e, e')L\{e') 

{e) + lLl^^{e)]<^r{^.^'mo-L\{e'))], 

(A.12) 



^-0 

oo 

1=0 



DP, 



where the "Legendre moment" P; := \ S-l Pi{l^'')f{p) 
of order I of the distribution function has been intro- 
duced. These Legendre moments can obviously be written 
as linear combinations of the angular moments Af,,,, :— 
^ J^j^ d/i which occur in the formulation of the 
transport equations used in our code. Detailed balance 
requires 



0r(e,6') = e-(-^')/^</.r(e,o. 



(A.13) 



fQj- Note that exploiting the in-out invariance of the 
transition rate (and therefore 0P*(e,e') ~ 0|"(e',e); 
e.g. [Cernohorsky 1994) leads to 



dee^B''°\e) ^Q, but 



dee^Bf\e)^0, 



(A.14) 

which means that neutrino number is conserved in the 
scattering process, whereas there is a nonvanishing en- 
ergy exchange between neutrinos and matter due to scat- 
terings. 



^ We assume spherical symmetry, whic h imp lies that the neu- 
trino distribution functions entering Eq. ( |A.6| ) are independent 
of the angle ip' . 
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Isoenergetic scattering: If for a particular scattering pro- 
cess the energy transfer between neutrinos and target par- 
ticles can be neglected, we have (/)°"'(e, e') = 0J"(e, e') =: 
(/);(e) • S{e — e'), and the source terms given by Eqs. (A. 10 



A.ll| , |A.12D simphfy to 

OO 



1=0 
,2, 



0, 

2^e2ii(e).(</,i(e)-0o(e)) 



(A.15) 
(A.16) 
(A.17) 



A. 1.3. Pair processes 



Applying the procedures outlined in Sect. |A.1.2| to the 
collision integral for the process of thermal emission 
and absorption of neutrino-antineutrino pairs by electron 
positron pairs and related reactions (e.g., nucleon-nucleon 
brcmsstrahlung) . the corresponding source terms read (see 
Brucnn|[l985|;|Pons et al.||l998|): 



/•OO 

Btp{€,h) = 2n de'e'^ 

OO 

{(1-/.) • ro{e,e')-Y,{2l + m{^^^n^^^W) 

OO 

+U ■ Y.{2l + l)Pl{^l)cj,'^{e,e')Ll{e'))], (A.18) 



i=0 



4°])(e) = 27r de'e''{</>P(e,e')-(l-io(e)-io(e')) 

OO 

^(2? + l)0f(e,e')iKe)^Ke')}, (A.19) 



(=0 



4^^(e) =2^ de' - 0P(e, e')^i(e) - 0?(e, e')^i(e) 

OO 

Y^{21 + l)0r(e,6')[(/ + l)ii+i(e) + lLi^i{e)]Li{e')] 

(A.20) 



1=0 



Here bars indicate quantities for antincutrinos. The ab- 
sorption and production kernels are related by 

0r(e,6') = (l-e(^+^')/^).</)P(6,O (A.21) 

in accordance with the requirement of detailed balance. 

A.2. Numerical implementation of various interaction 
processes 

In the following we describe the implementation of the 
various neutrino-matter interactions into our transport 
scheme. The expressions are mainly taken from the works 



of|Tubbs fc Schramnj (|1975D, [Schinder fc Shapiro| (|1982D, 



Brucnn (1985), and Mezzacappa & Brucnn (19931:,^ 

All average values of the source terms within individ- 
ual energy bins (cf. Eq. ^0|) are approximated to zeroth 
order by assuming the integrand to be a piecewise con- 
stant function of the neutrino energy e, and hence: 



B 



(fc) 

J + l/2 



BW(, 



J + 1/2 



(A.22) 



With the rest-mass energy of the electron, moC^ = 
0.511 MeV, and the characteristic cross section of weak in- 
teractions CTo := 4,{m^c^Gvf/{TT{hc)^) = 1.761-10-'''* cm^ 
(Gf is Fermi's constant), we define: 



(A.23) 



A. 2.1. Neutrino absorption, emission and scattering by 
free nucleons 

In the stand ard description of n e utrin o- nuclcon inter- 
actions (see Tubbs & Schra 





1975 




Brucnn 




1985 



Mezzacappa fc Brucnnj |1993q ), many-body effects for the 
nucleons in the dense medium, energy transfer between 
leptons and nucleons as well as nucleon thermal mo- 
tions are ignored. The corresponding simplifications al- 
low for a straightforward implementation of the processes, 
using Eq. ([A.5|) for the charged-current reactions and 
Eqs. ( A.15 - [AT7 ) for the neutral-current scatterings, re- 
spectively. Expressions for k* and the Legendre coeffi- 
cie nts 00(e ), 0i(£ ) can be computed from the rates given 
by Brucnn (1985) and Mezzacappa & Brucnn ( 1993c| ). 

For absorption of electron neutrinos by free neutrons 
{vq + n ^ c^ + p) the final result for the opacity is: 



2 , o 2^ 1--Fe-(£ + Q) 
1 - fu^{^) 



■{e + Q)^{e + QY~mlc^. (A.24) 

The opacity of the absorption of electron antincutrinos by 
free protons (f'c + P ^ c+ + n) is 



3.9i; 



l-fe+(£-Q ) 



Vpn 



(e-Q)V(^^ 



(A.25) 



if the neutrino energy e > rricC^ + Q, and K*(e) = 0, else. 
Here, Q :— m^c^ — m-p(? denotes the difference of the rest- 
mass energies of the neutron and the proton, and F^^ are 
the Fermi distribution functions of electrons or positrons. 
The constants — ^ and — 1-254 are the nucleon 
form factors for the vector and axial vector currents, re- 
spectively. The quantities rynp and 77p„ take into account 
the effects of fermion blocking of the nucleons. Making 
the approximations mentioned above, in particular ignor- 
ing the recoil of the nucleon, Brucnn (|198|, Eq. C14; see 
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also Reddy et al. 1998|) derived 



7np 



(27rfec)3^" 



{vc + A' ^ A + e^) is implemented by using Eq. (A. 5). 
For calculating the opacity of this process we adopt the 
description employed by Bruenn ( 1985| ) and Mezzacappa 



fc Brucnn (1993c) 



exp [Tpp -Tpn]-1 



(A.26) 



where rin, rip are the number densities of neutrons and 
protons, respectively. The corresponding degeneracy pa- 
rameters ipn and '0p are calculated by inverting the rela- 
tion (N stands for n or p) 



1 / 2 TONC^ T 
2^ \ [hcf 



3/2 



^1/2 (V'n) , 



(A.27) 



where F1/2 is the F ermi integral of index 1/2 (see, e.g., 
Lattimer & Swesty 1991, Sect. 2.2). The quantity rjp^ 
in Eq. ( A.25| ) is defined accordingly by exchanging the 
subscripts n and p in Eq. ( A.26| ). In the nondegenerate 
regime, where blocking is unimportant, one verifies the 
familiar limits rypn = rip and ?7np = n-a- 

Scattering of neutrinos off free nucleons (n or p) is me- 
diated by neutral currents only, which makes the distinc- 
tion between neutrinos and antineutrinos unneccessary. 
Neglecting nucleon recoil and nucleon thermal motions, 
the isoenergetic kernel obeys i?is(e,a;) cx (1 -t- w) (e.g. 
Bruenn 1985; Mezzacappa fc Bruenrj 1993c[ ), which im- 



plies (/);>i(e) = 0. The non-vanishing Legendre coefficients 
read 



00 (e) 



Btt' 



'?nn[l + igl] 

47ypp[(C7v-l)2 + fgi] 



01 (e) = — 



The factors 



Q J?7nn[l-5A] 



4,7pp[(C7v - If 



hi] 



(for n) , 
(for p) , 

(for n) , 
(for p) . 



'7nn/pp 



:= T 



an„/p 



(A.28) 
(A.29) 

(A.30) 



^Mn/p 

approximately account for nucleon final state blocking 
( Bruenn |1985D , and Cy = 1/2 2 sin^ with sin^ = 
0.23 for the Weinberg angle. In order to accurately re- 
produce the limits ?7nn/pp = "-n/p for nondegenerate (and 
nonrelativistic), and r^nn/pp = "n/p ' 37'/(2£^^/p) for de- 
generate (and nonrelativistic) nucleons (^n/p Fermi 
energy) we prefer using the analytic interpolation formula 
proposed by Mezzacappa fc Bruenn| ( |l993c) ) 



^iin/pp ^^n/p 



with 



Cn/p 



3r 

n/p 



instead of a direct numerical evaluation of Eq. ( |A.30|) 



iA.31) 



A. 2. 2. Neutrino absorption and emission by heavy 
nuclei 

As in the case of charged-current reactions with free nu- 
cleons, neutrino absorption and emission by heavy nuclei 



<{e)^g-gl-Np{Z)N^{N) 



Q') 



■{e + Q')^{e + Q'Y~mlc\ (A.32) 

if e > moC^ — Q', and K*(e) — else. Here Q' :— 
mA'C^ — mAC^ ~ MnC^ — AfpC^ + A denotes the energy dif- 
ference between the excited state A' = (N +1, Z — 1) and 
the ground state A = {N, Z). Following Bruenn ( 1985 ), we 
set A = 3 MeV for all nuclei. The internal structure of the 
nucleus with charge Z and neutron number N = A — Z is 
represented by the term 2/7 Np{Z)N\i{N), with the func- 
tions Np{Z), N\i{N) as given by Mezzacappa & Bruenn 
( |1993c) , Eqs. 31, 32). 

A. 2. 3. Coherent scattering of neutrinos off nuclei 

We use reaction rates for coherent neutrino scatterings 
off nuclei which include corrections due to the "nuclear 



1993c; Bruenn & Mezzacappa 



ion correlations (as described in Horowitz 1997). For a de 



1997D , and iorr 



form factor" (following an approximation by Mezzacappa 



tailed discussion of the numerical handling of both correc- 
tions, see the appendix of Bruenn & Mezzacappa (1997). 
The result of the latter reference can readily be used to 
calculate 5ig\ons' contribution of coherent scatter- 
ing to the rhs. of the first order moment equation. For 
the Boltzmann equation, we have to make a minor addi- 
tional approximation: The correction factor (^(e));^,^ as 
provided by Horowitz (1997) applies for the transport 

61 (which is 



opacity and thus for the combination 



^0 



given by Bruenn & Mezzacappa 1997) but not necessarily 
for the Legendre coefficients (f>Q and (pi individually, which 
are needed for calculating the collision integral i?is,ions 
for the Boltzmann equation. We nevertheless assume here 
that (pQ/iie) cx (^(e));^^^. This implies 



00/1 (e) -Ye^S-A' UAZiA, Z) ■ To/i{y) (5(e)), 



(A.33) 



with 

My) 
•^i(y) 
y 



-2y 



2y-l + e-^f 

9 ' 

y 

2-3y + 2y^ - {2 + y)e 

y3 

2(L07A)2/3e2(io-i3cni)2 
5 



(A.34) 



(he) 



containing the effects of the nuclear form factor, and 



ZiA,Z) 



[Ca - Cv) + (2 - Ca - Cv) 



2Z -A 



A 



(A.35) 
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Cv is defined as before (see Eqs. A.28 , A. 29 ) and Ca — 
1/2. In effect, the angular dependence of the ion- ion cor- 
relation correction is not exactly accounted for in the 
Boltzmann equation. The approximation, however, only 
influences the closure relation for the moment equations 
at a higher level of accuracy. It neither affects the Legendre 
coefficients at larger neutrino energies (where {S{e))^^^ is 
close to unity) nor does it destroy the consistency of the 
Boltzmann equation and its moment equations. Moreover, 
it has been shown that the overall effect of ion-ion corre- 
lations in stellar core collapse is only marginally relevant 



(Brucnn & Mezzacappa 1997; Ramp]: 2000) 



A. 2. 4. Inelastic scattering of neutrinos off charged 
leptons 

For the scattering of neutrinos off charged leptons, we ap- 
proximate the dependence of the scattering kernel on the 
scattering angle by truncating the Legendre expansions 
(Eqs. [A.IC , A. 11 , A. 12 ) at /„iax = 1- For our purposes, 



this has been shown to provide a sufficiently accurate ap- 



proximation (see 


3mit & Cernohorsky 


1996 




Mezzacappa 


fe Brucnn 


1993b 


), in particular also because the total 



scattering rate is correctly represented. Expressions for 
the Legendre coefficients are taken from the works of 



Yueh & Buchlei 


( 


1977); 


Mezzacappa & Brucnn 


Cernohorsky (1994 


); 


Smit & Cernohorsky (1996) 



He,e')=«iA}(e,e' 
with k e {I, II} and 



(A.36) 



[2'KhcY 7max(0,e'-e) 

(1 -F,{E, + e- e'))Hf{e, e', E,) . (A.37) 



The constants ai/n, which depend on the neutrino 
species and also the type of charged lepton (i.e. elec- 
tron or positron) involved in the scattering process, are 
combinations of the weak coupling constants (see, e.g.. 



Cernohorsky 1994). For the /i and t neutrinos and an- 
tineutrinos, which are treated as a single kind of neutrino 
in our code, we take the arithmetic mean of the corre- 
sponding coupling constants. For ^ = 0, 1 the functions 
gy"(e,e',£;c) (in units of MeV^) are given by |Yueh fc 



Buchler (1977, with the corrections noted by Brucnn 
1985|). 



Calculating the Legendre coefficients (e^ , e^/ ) is 

a computationally expensive task, since for all combina- 
tions (ej, Cji) and all radial grid points num erical integrals 
over the energy of the charged leptons (Eq. A.37 ) have to 
be carried out. It is, however, sufficient to explicitly com- 

n/out 



/ = /°'^) can be verified for the employed approximation 
and discretization of the collision integral and its angu- 
lar moments to within the roundoff error of the machine, 
provided that all integrals over neutrino energies are ap- 
proximated by simple zeroth-order quadrature formulae 
(cf. Eq. A. 22). Similarly, one can also verify the conserva- 
tion of particle number (Eq. A. 14) in the corresponding 
finite difference representation. 

In our neutrino transport code the dependence of the 
source terms on the neutrino distribution and its angu- 
lar moments is treated fully implicitly in time, i. e. th e 
time index of the corresponding quantities in Eqs. (|A.1C| - 
A.12| ) reads /""'"^ and i""*"^- The thermodynamic quanti- 
ties and the composition of the stellar matter which are 
needed to calculate the Legendre coefficients for inelas- 
tic scattering of neutrinos, however, are computed from 
£* and l^e*, representing the specific internal energy and 
electron fraction at the time level after the hydrodynamic 



substeps (cf. Sect. |3.6.lD . This allows one to save computer 
time since the Legendre coefficients must be computed 
only once at the beginning of each transport time step. 
Comparing to results obtained with a completely time- 
implicit implementation of the source terms (i.e. using 
e"+^ and l^"'"''^ for determining the stellar conditions that 
enter the computation of the Legendre coefficients) we 
have found no differences in the solutions during the core- 
collapse phase, where neutrino-electron scattering plays an 
important role in redistributing neutrinos in energy space 
and thus couples the neutrino transport to the evolution 
of the stellar fluid. 

A. 2. 5. Neutrino pair production/annihilation by 
thermal e"*" e~ pairs 



Following the suggestion of p^ons et al. ( 1998| ) we truncate 
the Legendre series (A.18-A.20) for the pair-production 
kernels at ^max — 2. The Legendre coefficients are given 
by ( |Bru( 



1985| ; pons et alj |1998D 

y2 



1 



1 



(27rfic) 



3 K*i(ej,e/) + a2*/(ej',ej)) : 

(A.38) 



with X := (e + e')/r and Q defined by Eq. ( |A.23D . Explicit 
expressions for the (dimcnsionless) functions {tj , eji ) 
and details about their efficient calculation are given by 
Pons et al. ( 1998| ). Also the coefficients ai, q;2 for the dif- 



ferent neutrino flavours can be found there (Pons et al 
1998| , Table 1) 



pute 0, ' [e,^) for e, < ey . The missing coefficients 2.6. Nucleon-nucleon bremsstrahlung 



can be obtained by exploiting a number of symmetry prop 



erties (see Cernohorsky 1994). In doing so, not only the 
computational work is reduced by almost a factor of two, 
but even more importantly, detailed balance (Bnes = 0, if 



Legendre coefficients for the nucleon-nucleon 
bremsstrahlung process, assuming nonrclativistic nucle- 
|ons, are calculated by using the rates given in Hannestad 
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& RafTelt (1998). The non- vanishing coefficients read: 



SCI 



(e---l) 



(A.39) 



with 



E 

{N} 
1 

^9 



0.25(r,g)2 



s{x) , 



(A.40) 



3/2 



T2 



37r2 



ttinc 



2mNc2 T 



Coulomb lattice corrections for the pressure, energy den- 
sity, entropy, and adiabatic index. 

In general, the regime of NSE is bounded from below 
by a complicated hyperplane in (p, T, y^)-space. For our 
purposes, however, it is sufficient to assume that the tran- 
sition between the two regimes (and thus the two equa- 
tions of state) takes place at a fixed value of the density 
po: If the density of a fluid element fulfils the condition 
p > po = 6 X 10^ gcm^'^, the temperature in a supernova 
core is usually large enough (T > 0.5 MeV) to ensure 
that NSE holds. Therefore the EoS of iLattimer & Swesty 



(A.41) 

In Eq. (A.41), pp — h - (Stt^ ^n)^^'^ denotes the Fermi mo- 
mentum for the nucleons, = 15 is the pion-nucleon 
"fine-structure constant" (see Hannestad & Raffclt 1998D , 
and Ton is the nucleon mass. The analytic form of the 
dimensionless fit functions g and s{x), which both de- 
pend on the nucleon density tin and temperature T of 
t he st ellar medium can b e foun d in planncstad fc Raffclt 
( |l998| ). The sum in Eq. (|A.39| ) runs over the individual 
proces ses vDim ^ nn; (N =^ n). vv'oi) ^ pp; (N — p) 



(1991) can be used at such conditions. If the density drops 



below the value pq, the non-NSE equation of state is in- 
voked. Since both EoSs employ a similar description of 
the physical properties of stellar matter at densities in the 
range between 10'' gcm~^and 10^ gcm~'^, the merging of 
them across the (T, ye)-plane at po is sufficiently smooth 
as far as, e.g., the pressure, internal energy density and 
chemical potentials as functions of density are concerned. 

A particular complication, however, arises with the 



chemical composition of the stellar plasma. Lattimcr & 



I ' ' ' 1 — I 1 — I — ' ' I ' 1 p 1 ■ ■ 

and t/i jnp ^ np: (N = np). The latter process is approx- Swestyj ( |199lD describe the baryonic part of the EoS as a 



imately taken into account by using the particle density 

n,, 



np ^/Ti^n-p and multiplying the corresponding contri- 
bution to the sum in Eq. (|A.3g| ) by a factor of 28/3 (cf. 
Thompson et al.|^000|). 



The production processes of neutrino-antineutrino 
pairs both by the annihilation of e"*" e~ pairs and by the 
nucleon-nucleon bremsstrahlung are implemented into the 
discretized equations (Eqs. ^ |2^, |2^) by applying the 
techniques described for the inelastic scattering reactions 
(see Sect. |A.2.4|). 



mixture of nucleons, a-particles and a heavy nucleus with 
in general non- integer mass and charge numbers {A,Z). 
The latter is considered to be representative of a mixture 
of heavy nuclei that coexist at given (p, T, 1^). Assuming 
NSE, the corresponding number densities of the nuclear 
constituents and the mass and charge numbers of the rep- 
resentative heavy nucleus are given as a function of the 
local values of p, T, and Y^. At the transition from NSE 
to non-NSE, the nuclear composition freezes out and one 
would ideally want to retain the baryonic abundances as 



Appendix B: Equation of state and nuclear 
burning 

In this section we give basic information about the equa- 
tion of state that is currently used in our simulations of 
stellar core collapse and supernovae. 

B.l. Numerical handling of the equation of state: 



We employ the equation of state (EoS) of Lattimer fc 



Swestyl ( |I991D with a value of Ks = 180 MeV for the 
incompressibility modulus of bulk nuclear matter. The 



other parameters can be found in the paper of Lattimer fc 



Swesty (1991). Since this EoS assumes matter to be in nu- 



clear statistical equilibrium (NSE) it is applicable only for 
sufficiently high temperatures and densities. A more gen- 
eral equation of state, which consistently extends to lower 
temperatures and densities not being available, we have 



supplemented the EoS of Lattimer fc Swesty (1991) with 
an EoS that describes a mixture of ideal gases of nucle- 
ons and nuclei with given abundances, plus ideal gases of 
electrons and positrons of arbitrary degeneracy and rela- 



tivity, plus photons (Janka 199£). The latter EoS includes 



given by the EoS of Lattimer & Swesty (1991) at po- While 
this is no problem in case of a Lagrangian hydrodynam- 
ics code, where the grid cells follow the evolution of indi- 
vidual fluid elements with specific information about the 
chemical composition attributed, the use of an Eulerian or 
moving grid requires a more complicated numerical pro- 
cedure. Here, additional advection equations (similar to 
Eq. ^ for the different nuclear components must be in- 
tegrated to trace the temporal evolution of the compo- 
sition on the whole numerical grid. With a finite, fixed 
number of such conservation equations to be solved in 
the code, one cannot allow for an arbitrarily large en- 
semble of different nuclear species with non-integer mass 
and charge numbers, as provided at freeze-out from NSE 
at po by the EoS of Lattimer fc Swesty ( 1991 ). Instead, 
we predefine a discrete set of nuclei {{Ak, Zk)}k=i,....N„^c 
at the beginning of the simulation, with k ~ 1 for neu- 
trons, k = 2 for protons, fc = 3 for a-particles and 
fc = 4, . . . , iViiuc for a number of suitably chosen heavy nu- 
clei. When the conditions in a computational cell change 
from NSE to non-NSE, we must map the NSE compo- 
sition nn(p o), np(po), na(po), n ( A(pn), z(p„)){po) as given by 
the EoS of Lattimer fc Swesty (1991) to our discrete sam- 
ple of species {{Ak, •^fc)}fe=i....,Af„uc with the constraints of 
charge neutrality and baryon number conservation. 
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The following procedure turned out to yield very satis- 
factory results. It is applied at each time step and in cells 
of the computational grid, which are close to a possible 
breakdown of NSE (i.e. which are close to po). Applying 
this procedure also in grid cells where NSE still holds and 
the EoS of Lattimer & Swesty (1991) is used, makes sure 
that the compositional information is available for a tran- 
sition to the non-NSE regime, if the freeze-out condition 
(p < po) is reached during a hydrodynamic time step. First 
we identify the densities of neutrons and a-particles with 
their NSE values: rtn = '^n(p), n-a — na{p). For given den- 
sity p and electron fraction Yc, the number densities of pro- 
tons and nuclei, rip and respectively, can then be 
determined from the requirements of charge neutrality and 
baryon number conservation. The index {Aj, Zj) points to 
a particular nucleus (Aj, Zj) e {{Ak, ^fe)}fe=4...Ar„„o, which 
is chosen from the set of non-NSE nuclei. It is associated 
with the representative heavy nucleus (A(p), Z{p)) accord- 
ing to the conditions: 



\Aj ^A\< \Ak -A\ yke{4,..., TVnuc} . 



Aj - Zj > A - Z , 
Zj<Z. 



(B.l) 



The first condition ensures that {Aj, Zj) is selected as the 
nucleus of the ensemble {Aj,Zj) S {{Ak, Zk)}k=i...N^^„ 
which is closest to {A{p),Z{p)) concerning its mass num- 
ber. A unique solution requires that no isobars exist in 
the ensemble. The second and third conditions guarantee 



B.2. Nuclear burning, dissociation and recombination 

Nuclear burning, dissociation of nuclei and the recombi- 
nation of free nucleons and a-particles to heavy nuclei can 
change the baryonic composition in the regime below the 
transition density (po = 6 x 10^ gcm~^) from the EoS of 
Lattimer fc Swesty] ( |1991| ) to the low-density EoS. We take 
into account complete burning of carbon, nitrogen, oxy- 
gen, magnesium and silicon as well as nuclear dissociation 
and recombination in approximate ways. 



Burning: If the temperature in a grid cell exceeds a thresh- 
old value characteristic of a certain nuclear reaction, the 
burning element is assumed to be instantaneously con- 
verted to the end product of the reaction. This is a reason- 
ably good approximation because the burning time scale, 
which is an extremely steep function of the temperature, 
is small compared with the time scale for hydrodynamic 
changes. Moreover, in dynamical supernova models one is 
mainly interested in the energy release associated with nu- 
clear reactions, but not in detailed information about the 
nucleosynthetic yields. 

When the condition for silicon burning applies, which 



is taken to be T ^ 4.5 x 10^ K ( |Hix fc Thielemani:| |1999 
Mezzacappa et al, 200l| ), we simply add up the local mass 



fractions of ^SSi and ^SNi: Xl^+^. = X^^^,, + XSa;, and 



56Ni -T ^>-28gi, 

0. Analogously, we treat the burning of ^^C 
to ^*Mg with a threshold temperature of 2.5 x 1 0^ K and 



set A^28gi 

24 



that 



rip > and n 



> 0. 



of 1^0, ^°N e, 24Mg to 28Si at T = 3.5 x 10^ K (|Woosley 
let al.||2002| ). 



Both the EoS of Lattimer fc Swesty (1991) and the 



non-baryonic part of our low-density equation of state are 
stored in tabular form for our calculations. Our table of 



the EoS of Lattimer fc Swesty (1991) has 180 and 120 



logarithmically spaced entries within the density range 
5 X 10^ gcm~'^ < p < 3 X 10^^ gcm^'^ and the tempera- 
ture range 0.1 MeV < T < 80 MeV, respectively, and 50 
equally spaced entries for 0.001 < < 0.6. Simil arly, the 
table for the lepton plus photon EoS ( Janka 1999) has 441 



entries for 10 gem < pYc < 10^^ gem and 141 en- 
tries for 8.6 X 10"^ MeV < T < 86 MeV. Intermediate val- 
ues are obtained by trilinear interpolation in log p, log T, 
Ye and bilinear interpolation in log{pYc), logT, respec- 
tively. 

Using EoSs in a discretized form we have performed 
the test calculation suggested by Swesty (199(: , §4) for ex- 
ploring the accuracy of the EoS evaluation. From our tests 
we can exclude any serious "unphysical entropy produc- 
tion or loss in otherwise adiabatic flows" (Swesty 1996| ). 
For several combinations of initial values for Y^ and s that 
are typical of conditions encountered in core-collapse sim- 
ulations, we found the deviations from adiabaticity to be 
negligibly small. The relative change of the entropy per 
baryon was less than lO^'^ for a density increase by a fac- 
tor of lO". 



Nuclear dissociation and recombination: When the density 
drops into the range p < po ^md the temperature fulfills 
T > 3x 10^ K at the same time the baryonic composition is 
adjusted to account for the photo-disintegration of nuclei 
and the recombination of free nucleons and a-particles to 
^^Ni. We distinguish between three different regimes (I- 
III), separated by curves pi{T),p2{T) in the p-T-plane 
(see Fig. pTh: 



log(pi(r)) := 11.62 + 1.5 log(T9) - 39.17/r9 , (B.2) 
log(p2(r)) := 10.60 + 1.5 log(T9) - 47.54/r9 , (B.3) 

where the logarithms are to base 10 and Tg := T/IO^K. 
The function pi (T) is derived from the equations of Saha 
equilibrium assuming a mixture of half (by mass) ^^Fe and 
half a-particles and free neutrons. Analogously, P2{T) is 
defined by half the mass being in a-particles and half in 
free nucleons (see Shapiro & Teukolsky 1983, Chap. 18). 



In detail the nuclear composition in the three regions is 
obtained from the following prescriptions: 

I) If p > pi(T), all free nucleons and a fraction < fi < 
1 of possibly existing a-particles in a corresponding 
grid cell are assumed to form ^^Ni. Pairing a suitable 
number of free neutrons and protons we therefore set 
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10- 



10 



10 



X 
X 



n+1 



max{X;'-X",0}, 



n+l 



with 



V* _ vn, f sr max{Afc - 2Zfc,0} 

— + /ll ■ 2^ ■ ^i- 



fc=4 



fe J 



fc=4 



fc J 



- V" _L f o niin{Afc - Zfc,Zfc} 



fe=4 



(B.5) 



The factor < /n < 1 is the degree of dissociation of 
heavy nuclei which is determined as outhned below. 
Ill) li p < p2 {T) , all heavy nuclei and a fraction < /m < 
1 of the a-particles is disintegrated into free nucleons: 



Fig. B.l. Area in the p-T- plane in the vicinity of the tran- 
sition between the EoS of Lattimer fc Swesty ( 1991 ) and 
the low-density EoS. The transition density po is indi- 
cated by the vertical dash-dotted line. The curves pi{T) 
and P2{T) bound three separate regions (HII) in the p- 
T-plane, where different prescriptions are used for deter- 
mining the nuclear composition in the low-density regime. 
Region I contains heavy nuclei and possibly free nucleons, 
in region II the composition is determined by a-particles 
and free nucleons, and in region III all nuclei and a- 
particles are dissolved. Below the temperature marked by 
the solid horizontal line the recombination of a-particles 
and free nucleons is suppressed. The threshold tempera- 
tures of the different burning reactions are given by the 
horizontal dotted lines. 



for the new composition: 



p -max{X;-X„",0}, 

^^Kc (1 - /l) • ^"rc ' 

= + /i • ^rHe + 2 • niin{X-, X;} . (B.4) 



III ■ 



A, 



+ /ill • -Xl\ 



He 



A-k — Zk 

k=4 

E^k 
'A'k 

fc=4 



^k 1 



^"ho - (1 ^ /ill) ■ ^"hc 1 

Vfc e {4, . . . , TVn 



(B.6) 



Nuclear transmutations release (or consume) nuclear bind- 
ing energy. Consequences of this effect should be taken 
into account in the hydrodynamical simulation. Since our 
EoS defines the internal energy density such that contri- 
butions from nuclear rest masses are included, the conver- 
sion between nuclear binding energy and thermal energy 
will automatically lead to corresponding changes of the 
temperature, pressure, entropy etc. of the stellar gas. 

The factors /i, /n and /ni are equal to unity in the 
corresponding regions I, II or III, respectively of the p-T- 
plane but vary between and 1 at the boundary curves 
(Eqs. B.2, B.3). This means that the composition changes 
take place only on the curves pi{T) and P2{T). The degree 
of dissociation or recombination (expressed by the values 
< /i,/ii,/iii < 1) is limited by the available internal 
energy. Starting out with given internal energy density 
(at a fixed value of p) and taking into account tempera- 
ture variations due to changes of the nuclear composition, 
the numerical algorithm maximizes /i, fu or /ni, while 
keeping the temperature on the boundary curves until the 
dissociation or recombination is complete. 



II) If p2(T) < p < Pi{T), a-particles are formed by the 
disintegration of heavy nuclei and the recombination 
of free nucleons. The new composition is given by: 



X"^ 



X 



1 =max{X* -X*,0}, 
i=max{X;-X„*,0}, 

x:^^ + 2-mm{x:,x;}, 



Xl+' = (1 - hi) ■ XJ! Vfc e {4, . . . , TVnuc} , 
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